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A  computerized  system  was  developed  to  control  an 
anaerobic  continuously  stirred  tank  reactor  utilizing 
cellulose  as  the  feedstock  for  different  operational 
objectives:  designated  methane  production  rate  (MPR)  and  the 
maximum  methane  production  rate  under  a  constraint  of 
minimal  chemical  oxygen  demand  (COD)  removal  efficiency. 

The  control  algorithm,  a  nonlinear  self-tuning 
regulator  (NSTR) ,  was  developed  to  direct  the  control 
action.   The  NSTR  was  the  combination  of  an  adaptive 
parameter  estimator  and  an  optimizer.   The  parameter 
estimator  used  a  complex  searching  technique  to  find  the 
best-fitted  parameters  based  on  historical  data  which  was 
updated  at  each  sampling  time.   With  the  parameters  found  in 
the  parameter  estimator,  the  optimizer  used  an  interval- 
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reducing  algorithm  to  locate  the  flow  rate  which  would 
produce  the  optimum  methane  production  rate  in  order  to 
achieve  the  defined  objective. 

To  verify  the  proposed  control  algorithm,  a  bench-scale 
reactor  was  set  up  and  equipped  with  hardware  including 
computer,  data  acquisition  and  control  system,  and  final 
control  elements.   The  parameter  estimator  and  optimizer 
were  coded  in  FORTRAN  and  executed  in  a  Sun  model  3/260 
minicomputer.   Based  on  the  information  obtained  from  the 
optimizer,  control  signals  were  generated  by  executing  the 
control  software  in  a  Zenith  model  248  microcomputer  and 
transferred  to  a  plug-in  power  relay  control  board  through  a 
Keithley  570  data  acquisition  and  control  workstation. 

In  one  of  the  designated  MPR  operations,  the  average 
controlled  output  was  0.4  L  CH4/L-day,  which  was  the  same  as 
the  desired  set  point.   The  results  for  another  set  point 
operation  showed  that  the  average  MPR  was  0.209  L  CH4/L-day, 
which  was  only  slightly  higher  than  the  desired  0.2  L  CH4/L- 
day.   For  the  maximum  MPR  operation,  the  results  showed  the 
reactor  performed  at  an  average  MPR  of  0.492  L  CH4/L-day  and 
remained  stable.   The  COD  removal  efficiency  in  all  three 
tests  was  higher  than  the  minimal  limit  of  85%. 
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CHAPTER  1 
INTRODUCTION 


Anaerobic  digestion  is  a  biological  process  for 
conversion  of  organic  matter  and  biomass  to  methane  and 
carbon  dioxide  in  the  absence  of  oxygen.   It  has  been 
recommended  and  implemented  for  waste  treatment  and  energy 
production  for  several  years.   This  process  involves  complex 
microbial  interactions  and  is  easily  harmed  by  many  factors. 
Therefore,  a  close  control  of  one  or  more  of  the  process 
variables  is  necessary  for  efficient  operation. 

Operation  and  control  of  conventional  anaerobic 
digesters  are  based  on  empirical  data  collected  with  the 
reactor  performing  at  steady  state.   These  data  are  not 
suitable  for  the  development  of  control  systems  that  handle 
dynamic  situations,  such  as  varying  feedstocks,  changing 
environmental  conditions,  or  even  modifying  operational 
objectives.   Additionally,  dynamic  control  requires  real 
time  data  of  variables  important  to  the  description  of  the 
anaerobic  digestion  process.   Some  of  the  variables  are 
difficult  or  too  costly  to  measure.   However,  a  computer  can 
be  used  to  estimate  the  unmeasurable  variables  through  the 
dynamic  simulation  of  the  process.   In  addition,  it  can 
conduct  data  acquisition  and  control  for  the  process. 
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As  computer  technology  develops  and  the  cost  becomes 
lower,  interest  in  computer-controlled  systems  is 
increasing.   Applications  of  computer  control  in  industries 
such  as  chemical  engineering  and  food  engineering  have  grown 
tremendously  in  recent  years.   However,  applications  in 
using  computer  control  in  anaerobic  digestion  processes  are 
very  limited,  and  they  have  not  involved  dynamic  control. 
Therefore,  development  of  a  computer-controlled  dynamic 
system  has  significant  meaning  at  this  stage. 

The  objectives  of  this  research  are  as  follows: 

1.  to  derive  a  dynamic  model  to  characterize  anaerobic 
reactor  performance  based  on  microbial  processes, 

2 .  to  develop  control  algorithms  for  operating  the  anaerobic 
reactor  at  a  designated  methane  production  rate  and  at 
the  maximum  methane  production  rate,  and 

3.  to  verify  the  proposed  control  algorithms  in  a  bench- 
scale  computer-controlled  system. 


CHAPTER  2 
LITERATURE  REVIEW 


Anaerobic  Process 


General  Background 

Methane  formation  from  organic  materials  has  been  known 
since  the  eighteenth  century  (Pine,  1971) .   However,  the 
application  of  this  process  was  not  reported  as  a  method  to 
treat  municipal  wastewater  until  1881  (McCarty,  1982) . 
Besides  its  ability  to  reduce  organic  pollution,  this 
process  also  produces  a  useful  by-product,  a  combustible 
gas,  methane.   As  a  result,  interest  in  this  technology  has 
grown  considerably  since  the  energy  crisis  occurred  in  the 
1970s.   This  process  can  be  applied  to  many  types  of  organic 
material,  including  animal  manures,  crop  residues,  food 
processing  wastes,  human  waste,  and  municipal  sludge,  for 
both  energy  production  and  waste  treatment. 

Process  Microbiology 

Anaerobic  digestion  is  a  microbial  process  which 
stabilizes  organic  matter  in  the  absence  of  oxygen  and 
produces  gases,  such  as  methane  and  carbon  dioxide,  and  a 
small  amount  of  sludge. 
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In  early  studies,  the  anaerobic  process  was  described 
as  a  two-stage  process  (McCarty,  1964a,  1964b) .   In  the 
first  stage,  complex  organics  were  hydrolyzed  and  fermented 
to  simple  organic  materials  by  a  group  of  facultative  and 
anaerobic  acid-forming  bacteria.   Most  of  the  end  products 
in  this  stage  were  organic  fatty  acids.   In  the  second 
stage,  a  group  of  complex  methane-forming  anaerobic  bacteria 
converted  the  products  from  the  acid-forming  stage  into  the 
gaseous  end  products,  methane  and  carbon  dioxide. 

Bryant  et  al.  (1967)  reported  that  a  nonmethanogenic 
species  produced  acetate  and  hydrogen  from  ethanol ,  and  the 
methanogenic  bacteria  only  degraded  methanol  and  acetate. 
Therefore,  another  group  of  bacteria  named  H2-producing 
acetogenic  bacteria  was  added  to  explain  this  process 

(Bryant,  1979) .   This  group  degraded  the  fatty  acids 
produced  from  the  first  stage  to  acetate,  C02,  and  H2. 
Recently,  another  group  called  homoacetogenic  bacteria  was 
found  to  synthesize  acetate  using  C02,  H2,  and  formate 

(Zeikus,  1980) .   However,  in  the  gastrointestinal  tract  of 
animals,  acetogenic  bacteria  were  probably  not  important 
because  of  the  short  retention  times,  and  only  the 
fermentative  bacteria  and  the  H2-utilizing  methanogenic 
bacteria  were  involved  in  the  partial  methane  fermentation 

(Hashimoto  et  al.,  1981;  Mclnerney  and  Bryant,  1981).   Flow 
diagram  of  the  microbial  transformations  in  the  anaerobic 
digestion  process  is  shown  in  Figure  2-1  (modified  from 
Zeikus,  1980) . 
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Figure  2-1, 


Flow  diagram  of  the  microbial  transformations 
in  the  anaerobic  digestion  process  (modified 
from  Zeikus,  1980) . 


Digester  Operation 

Since  the  anaerobic  process  involves  complex  biological 
reactions,  some  key  factors  which  should  be  considered  for 
successfully  operating  an  anaerobic  digester  include 
temperature,  retention  time  (both  hydraulic  and  solids) , 
mixing,  and  process  configuration  (Hashimoto  et  al.,  1981; 
Chynoweth  et  al.,  1984;  Loehr,  1984;  Parkin  and  Owen,  1986). 

Temperature.   Temperature  is  an  important  environmental 
factor  in  the  anaerobic  digestion  process  since  it  affects 
the  activity  of  the  microorganisms.   Methane  production  has 
been  found  at  temperatures  ranging  from  0  to  60  °C  (Kotze  et 
al.,  1969;  Svensson,  1984).   However,  practical  applications 
are  usually  performed  under  mesophilic  (20  to  40  °C)  or 
thermophilic  (50  to  60  °C)  conditions.   Within  each  range, 
the  maximum  methane  production  occurred  at  a  specific 
temperature  (Farguhar  and  Rovers,  1973) .   In  general,  higher 
temperature  will  promote  faster  reaction  rates  and  thus 
permit  operating  at  higher  loading  rates  without  reduction 
in  conversion  efficiency.   The  optimum  temperature  is  30  - 
37  °C  for  mesophilic  anaerobic  digestion  (Loehr,  1984)  and 
55  -  60  °C  for  thermophilic  condition  (Zinder  et  al.,  1984). 
Although  thermophilic  anaerobic  digestion  has  a  higher 
methane  production  rate  at  a  lower  hydraulic  retention  time 
(Shelef  et  al.,  1980;  Hill  et  al.,  1985;  Liao  and  Lo,  1985; 
Lo  et  al.,  1985),  thermophilic  conditions  are  not  widely 
used  because  of  the  high  requirement  for  external  energy  to 
maintain  the  desired  temperature  for  the  system.   Also,  the 


process  is  not  as  stable  as  under  mesophilic  conditions 
(Chynoweth  et  al.,  1984). 

Retention  time.   Retention  time  is  a  factor  which 
affects  bacterial  growth  and  the  system  cost.   It  can  be 
classified  as  hydraulic  retention  time  (HRT)  and  solids 
retention  time  (SRT) .   HRT  is  defined  as  the  time  required 
to  replace  the  fluid  volume  of  the  culture,  i.e.,  the 
reactor  volume  (V)  divided  by  the  flow  rate  (Q) .  When  flow 
rate  is  constant,  a  low  HRT  means  a  small  reactor  volume. 
Thus,  a  low  HRT  can  reduce  the  reactor  volume  as  well  as  the 
construction  cost  of  the  reactor.   Although  a  lower  HRT  is 
preferable,  the  washout  of  the  slow-growing  microorganisms 
may  cause  reactor  failure.   Therefore,  a  minimum  HRT  without 
sacrificing  the  process  stability  should  be  considered  in 
designing  a  reactor.   SRT  is  an  index  that  can  reflect  the 
detention  of  microorganisms  and  fixed  solids  in  the  reactor. 
Fundamentally,  SRT  is  defined  as  the  weight  of  solids  in  the 
system  divided  by  the  weight  of  solids  leaving  the  system 
per  unit  time  (Loehr,  1984).   If  SRT  is  less  than  the 
required  time  for  microbial  reproduction,  microorganisms 
will  be  removed  from  the  system  faster  than  they  can  be 
reproduced,  and  methane  formation  will  cease.   The  minimum 
SRT  for  anaerobic  digestion  was  estimated  in  the  range  of  2 
-  6  days  (Loehr,  1984) .   Many  researchers  reported  that  a 
longer  SRT  will  increase  the  stabilization  of  organic  matter 
and  enhance  the  stability  of  the  process  (McCarty,  1964b; 
Pfeffer,  1968;  Andrews,  1969).   Therefore,  increasing  the 
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SRT/HRT  ratio  is  desirable  when  designing  a  reactor  (Speece, 
1983)  . 

Mixing.   Mixing  is  used  to  (1)  maintain  an  intimate 
contact  between  the  microorganisms  and  their  substrates,  (2) 
efficiently  utilize  the  digester  volume,  (3)  disperse 
organics  and  inhibitory  substances  within  a  digester,  (4) 
prevent  stratification  and  temperature  gradients,  and  (5) 
minimize  scum  layer  formation  (Stafford,  1982;  Parkin  and 
Owen,  1986;  U.S.  Environmental  Protection  Agency,  1987).   In 
an  inefficiently  mixed  reactor,  zones  of  different  pH's  and 
different  temperatures  will  occur  and  reduce  the  efficiency 
of  the  reactor.   One  of  the  disadvantages  of  complete  mixing 
is  that  some  substrates  will  leave  the  reactor  unmetabolized 
(Kotze  et  al.,  1969).   The  quantity  lost  depends  on  the 
hydraulic  retention  time  and  the  conversion  efficiency  of 
that  substrate.   Fannin  et  al.  (1982)  showed  that  an  unmixed 
digester  had  a  higher  methane  production  rate  than  a 
continuously  stirred  tank  reactor  (CSTR)  using  sea  kelp  as 
substrate.   Mixing  should  be  considered  to  be  dependent  on 
the  type  of  anaerobic  process  and  the  type  of  substrate 
being  treated. 

Process  configuration.   Several  types  of  anaerobic 
processes  have  been  used  for  treating  food  processing 
wastes,  animal  manures,  wastewaters,  and  many  other  wastes. 
In  general,  they  can  be  classified  into  two  categories: 
single-stage  and  two-stage  processes.   In  a  single-stage 
process,  microbial  stabilization  and  liquid-solid  separation 
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occur  in  the  same  reactor.   Plug  flow  reactors  and  CSTR  are 
examples  of  this  type  of  process.   Some  modifications,  such 
as  baffles  for  the  plug  flow  system  and  mechanical  mixing 
for  the  CSTR,  were  added  to  increase  the  treatment 
efficiency  (McCarty,  1982). 

Since  it  was  realized  that  anaerobic  digestion  was  a 
two-stage  process  (acid  forming  and  methane  forming) ,  phase 
separation  was  thought  to  be  a  method  to  optimize  each  stage 
and  improve  the  efficiency  of  the  whole  system  (Smith  et 
al.,  1977;  Cohen  et  al.,  1979;  Ghosh  et  al.  1982;  Mata- 
Alvarez,  1987).   Many  studies  have  shown  that  two-phase 
anaerobic  digestion  can  operate  at  a  higher  organic  loading 
rate  and  have  a  comparable  or  higher  conversion  than  those 
in  single-stage  CSTR  digestion  (Colleran  et  al.,  1983;  Yang 
and  Chou,  1986;  Sarasevat  and  Khanna,  1986;  Ghosh,  1986; 
Stephenson  and  Lester,  1986;  Howerton  and  Young,  1987). 
Though  many  novel  and  promising  designs  have  been  developed, 
the  selection  of  an  anaerobic  process  depends  upon  the 
substrate  and  the  objectives  of  the  system. 

Process  Stability 

Anaerobic  digestion  processes  can  be  harmed  by  any 
factor  that  changes  the  environmental  conditions  or  causes 
the  accumulation  of  toxic  substances.   Examples  include 
organic  shock  loading,  abrupt  change  of  temperature,  and 
exposure  to  oxygen  or  antibiotics.   Some  key  variables  which 
can  reflect  the  stability  of  the  anaerobic  digestion  process 
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should  be  monitored,  controlled,  and  corrected  when  any 
adverse  condition  occurs.   These  variables  include  gas 
production  and  composition,  conversion  of  organic  matter, 
pH,  alkalinity,  and  volatile  acids  (Chynoweth  et  al.,  1984; 
Parkin  and  Owen,  1986) . 

Gas  production  and  composition.   The  gaseous  products 
of  anaerobic  digestion  are  methane,  carbon  dioxide,  and  some 
trace  gases  such  as  hydrogen  and  hydrogen  sulfide.   Although 
gas  production  is  a  direct  index  for  the  conversion  of 
organic  matter,  the  destruction  of  organic  matter  may 
contribute  to  the  formation  of  gases  other  than  methane. 
Therefore,  methane  production  is  a  more  sensitive  indicator 
for  overall  digester  performance.   Typically,  the 
destruction  of  1  gram  of  COD,  regardless  of  substrate 
source,  will  produce  0.35  std  liter  of  methane  (McCarty, 
1964a) .   When  methane  content  decreases  and  carbon  dioxide 
increases,  it  indicates  digester  instability  or  inhibition 
of  the  methane  bacteria. 

Conversion  of  organic  matter.   To  determine  the 
efficiency  of  the  digestion,  measurements  of  organic  matter 
are  required.   The  commonly  used  measurements  of  organic 
matter  include  volatile  solids  (VS)  and  chemical  oxygen 
demand  (COD) .   Other  measurements  such  as  5-day  biochemical 
oxygen  demand  (BOD5)  and  total  organic  carbon  (TOC)  are  also 
used  as  a  measure  of  organic  content.   Because  the  BOD5  test 
requires  a  long  time  to  get  the  result,  it  is  not  suitable 
for  reactor  control.   Although  TOC  is  a  more  rapid 
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measurement,  the  cost  of  the  instrumentation  is  high  when 
compared  with  VS  and  COD.   If  the  organic  matter  in  the 
influent  and  effluent  are  measured,  the  removal  efficiency 
can  be  determined.   When  the  reactor  is  in  steady  state,  the 
removal  efficiency  of  organic  matter  will  be  stable.   In 
addition,  measurements  of  the  influent  and  effluent  organic 
levels  are  necessary  for  calculation  of  the  system  mass 
balance. 

Alkalinity,  volatile  acids,  and  pH.   Alkalinity, 
volatile  acids  (VA) ,  and  pH  are  interdependent  and  very 
important  in  control  of  a  reactor  (McCarty,  1964c;  Kroeker 
et  al.,  1979;  Asinari  Di  San  Marzano  et  al.,  1981;  Stafford, 
1982;  Parkin  and  Owen,  1986).   The  optimum  pH  for  anaerobic 
digestion  is  between  6.6  and  7.6  (McCarty,  1964c).   When  pH 
deviates  significantly  from  this  range,  it  indicates  the 
imbalance  or  failure  of  the  digester.   Alkalinity  is  the 
buffering  capacity  of  the  digester.   Under  stable 
conditions,  bicarbonate  alkalinity  is  approximately  equal  to 
total  alkalinity.   When  volatile  acids  begin  to  increase, 
bicarbonate  alkalinity  is  neutralized,  and  volatile  acid 
alkalinity  results.   As  the  bicarbonate  alkalinity  is 
depleted  by  the  accumulation  of  the  volatile  acids,  the  pH 
decreases  and  the  condition  of  the  digester  will  become 
toxic  to  methane  bacteria.   If  the  bicarbonate  alkalinity  is 
in  the  range  of  about  2500  -  5000  mg/L,  even  a  large 
increase  in  volatile  acids  can  be  handled  (McCarty,  1964c) . 
Volatile  acids  such  as  acetic,  propionic  and  butyric  acids 
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are  formed  by  acid-forming  bacteria.   The  most  important 
volatile  acid  in  methane  fermentation  is  acetic  acid,  which 
contributes  72%  of  methane  formation  (McCarty,  1964b) .   In  a 
healthy  reactor,  the  formation  of  volatile  acids  and  the 
utilization  by  the  methane-forming  bacteria  is  balanced. 
Any  factor  that  stimulates  the  production  of  acids  or 
inhibits  the  methane- forming  bacteria  will  cause  the 
accumulation  of  the  volatile  acids,  and  thus  decrease  the  pH 
and  alkalinity. 

Kinetics  and  Modeling 

Kinetics  and  modeling  are  essential  for  process  design, 
performance  prediction,  operation,  and  control.   With  the 
help  of  the  digital  computer,  hypotheses  of  the  biological 
process  can  be  tested  more  economically  and  faster  than 
before.   Additionally,  the  feedback  of  the  calculated 
results  can  help  scientists  to  prove  or  modify  their 
understanding  of  the  experimental  data,  and  another 
experiment  can  be  planned  and  initiated.   The  components  in 
the  loop  of  experimentation,  observation,  hypothesis,  and 
modeling  are  interdependent  and  can  be  beneficial  to  each 
other. 

General  Kinetic  Model 

Aiba  et  al.  (1973)  described  kinetics  as  the  rates  of 
cell  synthesis  and/or  fermentation  product  formation  and  the 
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effect  of  environment  on  these  rates.  The  kinetic  equations 
of  bacterial  growth  were  introduced  by  Monod  (1942) .  Later, 
researchers  developed  some  modifications  of  "Monod  kinetics" 
(Lawrence  and  McCarty,  1969;  Contois,  1959;  Andrews,  1968). 

Monod  kinetics.   In  1942,  Monod  proposed  a  relationship 
between  the  limiting  substrate  concentration  and  the 
microbial  growth  rate.   Since  then,  it  has  been  used  as  the 
basis  for  most  of  the  anaerobic  digestion  models.   With  the 
same  form  as  the  rate  equation  for  a  one-substrate  enzyme- 
catalyzed  reaction  (Lehninger,  1982) ,  the  Monod  equation 
states  that 


H   =   (2-1) 

Ks  +  S 

where     /x  =  specific  growth  rate,  time"1 

Mm  =  maximum  specific  growth  rate,  time" 

S  =  substrate  concentration,  mass/volume 

Ks  =  saturation  constant,  S  at  jli  =  M„/2. 


Substrate  utilization  rate  and  decay  coefficient.   The 
Monod  equation  for  the  substrate  utilization  rate  is 


dF     k'M'S 
=  (2-2) 


dt      Ks  +  S 

where     dF/dt  =  substrate  utilization  rate, 

mass/volume-time 
M     =  microorganism  concentration,  mass/volume 
k     =  maximum  specific  rate  of  substrate 

utilization,  time"1. 

Including  a  decay  rate  (Lawrence  and  McCarty,  1969) ,  the  net 
growth  rate  can  be  expressed  as 
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dM         dF 

=  a* ( )  -  b«M  (2-3) 


dt        dt 


where     dM/dt  =  microorganism  net  growth  rate, 

mass/volume-time 
a     =  growth  yield  coefficient,  dimensionless 
b     =  microorganism  decay  coefficient,  time  . 


Decrease  in  cell  mass  may  be  caused  by  the  energy  required 
for  cell  maintenance,  death  and  predation  of  the  cells. 
Therefore,  the  decay  coefficient  b  in  the  above  expression 
is  a  lumped  factor  (Metcalf  &  Eddy,  1979) . 

Combining  the  above  two  equations,  a  relationship 
between  specific  growth  rate  and  specific  substrate 
utilization  rate  was  obtained: 


dM/dt      a«k«S 
= b  (2-4) 


M         Ks  +  S 

The  quantity  (dM/dt) /M  is  the  specific  growth  rate  ju,  which 
is  equal  to  the  mean  hydraulic  retention  time  under  steady 
state  condition  for  a  completely  mixed  reactor.   And  a«k  is 
equal  to  the  maximum  specific  growth  rate  /im. 

Contois  kinetics.   Contois  (1959)  presented  a  model 
which  related  the  bacteria  specific  growth  rate  to  the 
microorganism  concentration  and  the  limiting  substrate 
concentration : 


M  =  (2-5) 

BX  +  S 
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where     X    =  microorganism  concentration,  mass/volume 
B    =  dimensionless  kinetic  parameter. 


In  his  report,  Contois  pointed  out  that  the  saturation 
constant  Ks  in  the  Monod  equation  is  a  function  of 
microorganism  concentration. 

Substrate  inhibition  kinetics.   In  order  to  present  the 
inhibitory  effect  of  a  high  concentration  of  substrate, 
Andrews  (1968)  proposed  a  substrate  inhibition  model: 


/x  =   (2-6) 

1  +  K  /S  +  S/K, 


where  K.  =  inhibition  coefficient,  mass/volume. 

According  to  Andrews,  his  model  was  able  to  predict  process 
failure  due  to  organic  and  hydraulic  overloading. 

First-order  kinetics.   After  being  unsuccessful  in 
fitting  his  experimental  data  with  the  Monod  equation, 
Pfeffer  (1974)  adopted  the  first-order  rate  expression: 


dS 
=  -K-S  (2-7) 


dt 

where  K  =  rate  constant,  time"1  (not  the  same  as  Ks  or  k)  . 

Using  an  estimated  maximum  gas  yield  (0.547  liter  gas/g  VS 
added) ,  Pfeffer  calculated  the  substrate  concentration  and 
deduced  rate  constants. 

Other  forms  of  kinetic  models.   Some  other  expressions 
of  growth  kinetics  have  also  been  developed,  such  as  product 
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inhibition  kinetics  by  Aiba  et  al.  (1968),  multi-step  Monod 
equation  by  Hill  (1982c),  Hill  et  al.  (1987)  and  Braha  and 
Hafner  (1985) ,  and  some  exponent-included  curve  fitting 
models  (Bailey  and  Ollis,  1986) . 

Anaerobic  Digestion  Modeling 

Anaerobic  digestion  models  can  be  classified  into  two 
categories:  steady  state  models  and  dynamic  models.   A 
steady  state  model  is  useful  in  predicting  gas  production, 
gas  composition  and  effluent  characteristics  under  stable 
operating  conditions  of  feedstock,  temperature,  and  loading 
rate.   However,  a  steady  state  model  is  not  suitable  for 
predicting  failure  and  transient  behavior  of  the  process.   A 
dynamic  model  can  be  used  to  make  these  predictions.   Thus, 
a  dynamic  model  can  be  useful  in  predicting  response 
following  start-up  and  shock  loading. 

Andrews  (1968,  1969)  showed  that  a  pure  Monod  equation 
is  not  valid  in  a  reaction  in  which  volatile  acids  serve  as 
the  substrate  and  in  which  they  are  also  inhibitory  to  the 
methane  production  when  they  accumulate  to  a  toxic  level. 
He  modified  the  Monod  equation  and  presented  a  substrate 
inhibition  model  (eq.  2-6).   In  this  modified  Monod  model, 
Andrews  used  an  inhibition  function  to  relate  volatile  acids 
concentration  and  specific  growth  rate  for  the  methane 
bacteria.   Also,  unionized  volatile  acids  were  considered  as 
both  the  limiting  substrate  and  the  inhibiting  agent. 
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Using  the  same  inhibition  kinetics,  Andrews  and  Graef 
(1971)  modified  their  previous  model  by  adding  the 
interactions  between  the  liquid,  gas,  and  biological  phase 
of  the  digester  to  predict  the  volatile  acids  concentration, 
alkalinity,  pH,  gas  flow  rate,  and  gas  composition,   the 
same  kinetic  model  was  adopted  by  Bolle  et  al.  (1986)  and 
Moletta  et  al.  (1986)  in  their  studies  on  different 
substrates. 

Hill  and  Barth  (1977)  extended  Andrews'  model  to 
include  the  inhibitory  effect  of  unionized  ammonia  on  the 
methanogens  and  developed  a  simulation  model  for  digestion 
of  animal  wastes.   Using  the  dimensionless  form  of  coupled 
time-dependent  eguations,  Van  den  Heuvel  and  Zoetemeyer 
(1982)  applied  substrate  inhibition  kinetics  to  a  CSTR  with 
cell  recycle  to  predict  the  process  behaviors  in  both  steady 
state  and  dynamic  state  and  to  establish  the  criteria  for 
reactor  operation. 

The  above  modified  Monod  kinetic  models  are  accurate  in 
predicting  process  failure  and  optima  but  have  difficulty  in 
identifying  the  problem  parameters.   Contrary  to  the 
complexity  of  the  Monod  type  of  model,  the  first-order 
kinetic  model  is  simple  in  terms  of  input  parameters  (Grady 
et  al.,  1972,  Pfeffer,  1974,  Oleszkiewicz  and  Koziarski, 
1986) .   However,  this  type  of  model  is  unable  to  predict  the 
optima  and  process  failure  (Hill,  1983). 

Another  type  of  model  was  developed  by  Chen  and 
Hashimoto  (1978,  1980)  by  modifying  the  Contois  model.   They 
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proposed  the  following  equation  to  relate  the  influent  and 
effluent  substrate  concentration  to  the  specific  growth  rate 
for  a  completely  mixed  reactor  without  solids  recycle: 


*Vs/so 

H    =   (2-8) 

(l-K)-S 
K  +  


where  S0  and  S  are  the  influent  and  effluent  substrate 
concentration  (mass/volume)  respectively,  and  K  (different 
from  the  K  in  the  first-order  kinetics)  is  a  dimensionless 
kinetic  parameter.   In  this  model,  the  input  parameters 
reduce  to  only  /im  and  K.   Incorporating  the  relationship 
between  the  destroyed  substrate  and  the  theoretical  methane 
production,  this  model  can  be  used  to  predict  effluent 
substrate  concentration,  methane  production  rate  and  the 
optimal  retention  time  for  the  maximum  methane  production 
rate,  and  the  maximum  substrate  utilization  rate.   To  obtain 
the  maximum  specific  growth  rate  nm   for  the  above  model, 
Hashimoto  et  al .  (1981)  used  a  linear  relationship  between 
temperature  and  nm: 

Hm   =    0.013(T)  -  0.129  (2-9) 

where  T  is  the  temperature  between  20  and  60  °C. 

Using  the  above  model,  the  effects  of  temperature, 
influent  volatile  solids  concentration  and  hydraulic 
retention  time  on  the  kinetic  constant  K  were  evaluated 
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(Hashimoto  and  Hruska,  1982;  Hashimoto  and  Hruska,  1984). 
It  was  found  that  K  increased  exponentially  as  influent 
substrate  concentration  S0  increased  and  was  described  by 

K  =  0.8  +  0.0016«exp(0.06'So)  (2-10) 

for  cattle  manure,  and 

K  =  0.6  +  0.0206«exp(0.51»So)  (2-11) 

for  swine  manure. 

Chen  and  Hashimoto's  model  has  been  adopted  in  many 
studies.   In  anaerobic  treatment  of  landfill  leachate,  Lema 
and  Ibanez  (1985)  found  kinetic  parameters:  nm  =   0.275  d*  ,  K 
=  0.4  65  and  the  refractory  portion  R  =  0.1.   Mata-Alvarez 
and  Llabres  (1988)  successfully  applied  the  model  to  down- 
flow  stationary  fixed  film  (DSFF)  reactors  for  swine  waste 
digestion  and  found  kinetic  parameters:  /xm  =  0.1  d"  and  K  = 
0.9.   Hill  (1982a)  developed  design  criteria  for  the  maximum 
methane  production  in  an  animal  waste  digestion  system. 
Hill  (1982b)  further  modified  the  model  by  subtracting  the 
death  coefficient  Kd  as  10%  of  the  maximum  specific  growth 
rate  (fim)    from  the  net  specific  growth  rate  fj.   and 
established  the  optimum  operational  criteria  for  the 
digestion  of  animal  manure. 

Although  the  above  modified  Contois  model  has  the 
advantages  of  simplified  parameter  input  and  can  predict 
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failure  due  to  washout  of  the  microorganisms,  it  is  only 
applicable  to  steady  state  conditions  and  unable  to  predict 
process  failure  due  to  inhibition  of  the  methanogens  (Hill, 
1983)  . 

In  an  attempt  to  achieve  both  accuracy  and  simplicity 
for  the  dynamic  process,  Hill  (1983)  developed  a  "lumped 
parameter"  method  based  on  Monod  kinetics  for  modeling  a 
continuous  expanding  reactor  (CER) .   The  uniqueness  of  this 
method  is  in  lumping  two  or  three  basic  Monod  parameters 
into  one  parameter  that  varies  with  waste  type.   In  this 
method,  any  type  of  waste  can  be  categorized  by  two  factors: 
biodegradability  and  acid  factor.   The  biodegradability  is 
defined  as  the  portion  that  the  volatile  solids  can  be 
destroyed  at  an  infinite  detention  time,  and  the  acid  factor 
is  the  fraction  of  the  biodegradable  volatile  solids  that 
are  in  the  form  of  volatile  acids.   Other  than  these  two 
factors,  the  other  parameters  used  in  the  model  are 
constant,  regardless  of  the  waste  type.   In  order  to  better 
predict  the  transient  behavior  of  the  process,  Hill  et  al. 
(1983)  modified  the  above  "lumped  parameter"  model  by 
changing  the  form  of  the  death  coefficient  Kd  from  a 
constant  (0.1  nm)    to  a  Monod-based  function. 

The  models  described  above  are  two-microorganism 
models,  including  only  acetogens  and  methanogens.   Hill 
(1982c)  developed  a  dynamic  model  with  four  groups  of 
microorganisms,  including  acetic  acid  metabolism 
(hydrogenogenic  bacteria)  and  the  reduction  of  carbon 
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dioxide  with  hydrogen  (homoacetogenic  bacteria) .   In  this 
model,  kinetic  constants  were  developed  from  computer 
iteration  and  basic  stoichiometry. 

Though  much  research  has  been  done,  the  current 
modeling  efforts  are  concentrating  on  the  quantifying  of  the 
unidentified  microorganisms.   With  better  understanding  of 
the  microorganisms  involved  and  their  reactions  in  the  basic 
process,  a  more  sophisticated  model  can  be  developed  which 
will  better  predict  process  performance. 

Computer-Controlled  Anaerobic  Digestion  Process 

Since  anaerobic  digestion  processes  are  easily  harmed 
by  a  change  in  environmental  factors  and  the  accumulation  of 
toxic  substances,  close  control  of  this  process  is  necessary 
to  keep  the  reactor  healthy  and  to  achieve  a  specific  goal 
such  as  maintaining  a  designated  methane  production  rate.   A 
digital  computer  can  be  implemented  to  achieve  the  above 
objectives  due  to  its  high-speed  computation  and  large 
information  storage  capacity.   In  addition,  the  cost  of 
microcomputers  has  decreased  dramatically  in  recent  years, 
and  this  makes  computer  control  more  attractive  than  ever. 

Computer  control  of  the  anaerobic  digestion  process  is 
a  field  which  combines  the  expertise  of  process 
biochemistry,  control  theory,  sensor  development  and 
computer  technology.   General  control  theory  has  been  well 
established  and  applied  in  conventional  chemical  and 
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electrical  engineering  processes  and  can  be  adapted  to  this 
particular  process.   Though  sensors  for  direct  measurements 
of  cell  concentration  and  activity  are  not  available,  this 
problem  can  be  solved  with  the  help  of  mathematical  modeling 
and  monitoring  of  indirect  measurements  (Wang,  1984;  Shimizu 
et  al.,  1984;  Goto,  1986).   Following  the  development  of 
computer  technology,  the  application  of  the  computer  in 
process  control  has  progressed  from  the  earlier  analog 
controller  to  direct  digital  control  (DDC) .   Also,  functions 
of  control  systems  have  improved  from  data  acguisition  to 
data  analyzing  and  decision  making. 

Computers  emerged  in  control  systems  for  missiles  and 
aircraft  around  1950.   Their  application  in  industrial 
process  control  has  been  increasing  exponentially  since  1960 
(Astrom  and  Wittenmark,  1984) .   Along  with  this  trend, 
computer  control  has  been  used  extensively  in  conventional 
activated  sludge  plants  (Horan  et  al.,  1985;  Crowther  et 
al.,  1977;  Corder  and  Lee,  1986;  Vasicek,  1982;  Arthur, 
1982)  and  aerobic  processes  (Van  Breugel  et  al.,  1986). 

Computer  control  has  been  applied  in  only  a  few 
anaerobic  digestion  process.   Guarino  (1972)  used  a 
radioactive-type  density  gauge  to  measure  the  percent  solids 
in  sludge  in  order  to  meet  the  desired  solids  level  for 
sludge  treatment.   In  addition,  C02  was  monitored  on-line 
and  pH,  volatile  acids,  and  alkalinity  were  determined  in 
the  laboratory  and  supplied  as  the  reference  indices  for 
control. 
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Graef  and  Andrews  (1972)  reported  a  control  strategy, 
scrubbing  of  carbon  dioxide  from  the  biogas  with  subseguent 
recycle  of  the  treated  gas,  that  could  provide  the 
adjustment  of  digester  pH  and  thus  remove  the  carbonic  acid 
without  adding  base.   Andrews  (1978)  also  proposed  another 
control  strategy  which  recycled  the  sludge  from  a  second 
stage  digester  and  used  methane  production  rate  as  a 
feedback  signal.   Though  the  control  action  was  not  new, 
methane  production  rate  as  the  indicator  of  the  digester 
condition  was  a  novel  idea.   Both  the  above  strategies  have 
been  proven  feasible  through  simulation  studies. 

Rozzi  et  al.  (198  3)  developed  a  method  for  monitoring 
the  carbon  dioxide  content  of  biogas  by  measuring  the  pH  of 
a  sodium  bicarbonate  solution  which  the  biogas  was  passed 
through.   This  method  had  to  be  calibrated  with  a  standard 
N2/C02  gas  mixture  to  arrive  at  a  temperature- independent 
condition. 

An  automatic  control  system  for  anaerobic  sludge 
treatment  was  developed  by  Russell  et  al.  (1985).    This 
system  was  designed  primarily  for  maintaining  steady  state 
of  the  process.   Variables  controlled  in  this  system 
included  temperature,  pH,  flow,  and  tank  level.   The  1800  m3 
anaerobic  digester  had  a  93%  removal  for  soluble  COD  and 
soluble  BOD,  and  produced  a  total  gas  flow  of  3.8  L/L-day 
with  78%  of  methane. 

Whitmore  et  al.  (1986,  1987)  demonstrated  a  control 
system  which  monitored  the  dissolved  hydrogen  through  a  mass 
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spectrometer.   The  hydrogen  signal  from  the  mass 
spectrometer  in  a  feedback  loop  was  used  to  regulate  the 
hydrogen  level  by  controlled  addition  of  the  carbon  source 
(glucose) .   The  results  showed  that  volumetric  methane 
production  rate  of  1.4  L  CH4/L-day  could  be  maintained  when 
hydrogen  concentration  was  controlled  at  1  jiM  for  a 
mesophilic  anaerobic  digester. 

All  systems  listed  above  were  for  constant  set-point 
control .   They  are  suitable  only  for  the  designated 
feedstock  at  steady  state  but  not  proper  for  dynamic 
control.   A  reliable,  efficient  anaerobic  digestion  process 
could  be  expected  if  an  integrated  system  of  dynamic 
modeling  and  computer  control  can  be  developed. 


CHAPTER  3 
CONTROL  SYSTEM  DEVELOPMENT 


The  major  components  in  a  typical  control  system 
include  the  process,  measuring  devices,  controller,  final 
control  element  and  transmission  lines  for  measurements,  and 
the  control  signal  (Figure  3-1) .   This  chapter  will  discuss 
the  development  of  each  component,  including  hardware  and 
software. 

Process 

Reactor  Setup 

A  continuously  stirred  tank  reactor  (CSTR)  with  10- 
liter  liquid  volume,  shown  in  Figure  3-2,  was  set  up  for 
this  study.   The  reactor  was  mixed  intermittently  for  20 
minutes  every  hour  with  a  Dayton  2Z814  motor  (Dayton 
Electric  Mfg.,  Co.,  Chicago,  Illinois)  at  a  speed  of  33  rpm. 
During  twice  per  day  feeding,  the  feedstock  was  kept  in  a 
glass  bottle  and  was  also  mixed  for  the  same  length  of  time 
as  the  reactor  by  using  a  magnetic  stirrer.   In  order  to 
eliminate  the  effect  of  temperature  variations  on  the 
process,  the  reactor  was  submerged  in  an  insulated  water 
bath  where  the  temperature  was  maintained  at  37  ±  1  °C  by 
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using  a  thermostat-controlled  copper  tubing  immersion 
heater. 

Influent  and  effluent  flow  were  controlled  by  using 
Masterflex  pumps  (Cole-Parmer  Instrument  Co.,  Chicago, 
Illinois)  at  a  rate  of  15  ml/min.   To  insure  that  the 
influent  and  effluent  were  well  mixed,  the  flow  pumps  were 
started  5  minutes  after  mixing  of  the  reactor  liquids  and 
influent  were  initiated. 

Gas  produced  from  the  reactor  was  passed  through  a 
water-trap  bottle  to  a  sampling  bulb  for  gas  composition 
measurement  and  through  a  gas  meter  for  the  measurement  of 
gas  volume. 

Feedstock  Preparation 

In  order  to  have  better  control  of  the  mass  flow,  a 
synthetic  substrate,  a  purified  and  partially  depolymerized 
cellulose  (Avicel,  FMC  Corporation,  Philadelphia, 
Pennsylvania) ,  was  chosen  as  the  primary  component  of  the 
feedstock.   The  inoculum  used  for  this  experiment  was 
effluent  from  an  anaerobic  digester  using  sorghum  as 
feedstock.   To  start  the  experiment,  9  liters  of  the 
inoculum  was  mixed  with  1  liter  of  the  feedstock  and  fed 
into  the  reactor. 

Nutrients  and  vitamins  added  to  the  feedstock  were 
modified  from  the  studies  of  Owen  et  al.  (1979)  and  were 
prepared  as  concentrated  stock  solutions.   The  prepared 
feedstock  and  stock  solutions  were  stored  at  4  °C.   The 
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stock  solutions  and  the  procedure  for  preparing  feedstock 
are  listed  in  Tables  3-1  and  3-2.  The  characteristics  of 
the  inoculum  and  the  feedstock  are  presented  in  Table  3-3. 


Table  3-1.   Stock  solutions  for  preparation  of 
feedstock. 


Solution 

Compound 

Concentration    (g/L) 

Sol-3 

(NH4)2HP04 

26.7 

Sol-4 

CaCl2«2H20 

16.7 

NH4C1 

26.6 

MgCl2«6H20 

120.0 

KC1 

86.7 

MnCl2-4H20 

1.33 

CoCl2«6H20 

2.0 

H3BO3 

0.38 

CuCl2«2H20 

0.18 

Na2Mo04«2H20 

0.17 

ZnCl2 

0.14 

NiCl2«6H20 

0.15 

H2W04 

0.007 

Table  3-2.   Procedure  for  preparing  feedstock  at  a 
concentration  of  30  g  cellulose/L*. 


Step 


Instruction 


1. 
2. 


3. 
4. 


Prepare  a  2-liter  flask. 
Add:  2.8  ml  of  Sol-3 

14   ml   of   Sol-4 

6.98    g   of   NaHC03 

3.33  g  of  yeast  extract 

3.33  g  of  casein 

30  g  of  cellulose  (Avicel) 
Add  distilled  water  up  to 

1  liter. 
Store  at  4  °C  until  use. 


*  To  prepare  22.5  g/L  cellulose  feedstock,  three 
portions  of  30  g/L  cellulose  solution  was  mixed 
with  one  portion  of  distilled  water. 
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Table  3-3.   Characteristics  of  the  inoculum  and 
feedstock. 


Characteristics 

Inoculum 

Feedstock 

PH 

7.16 

7.75 

TS,  g/L 

20.0 

37.2 

VS,  g/L 

12.4 

30.8 

COD,  g/L 

13.3 

40.3 

VA,  mg/L 

as 

130 

510 

acetic 

acid 

*  Average  value  for  16  different  lots  of  feedstock 
at  a  concentration  of  30  g/L  cellulose.   For  finding 
the  characteristics  of  22.5  g/L  cellulose  feedstock, 
a  conversion  factor  of  0.75  need  be  multiplied. 

Operation 

According  to  the  mode  of  operation,  there  were  two 
phases  for  this  study:  the  first  phase  included  start-up  and 
steady  state  operation;  the  second  phase  was  the  dynamic 
state  operation  for  verifying  the  control  algorithm.   To 
start  the  reactor,  a  feedstock  with  a  concentration  of  3  0 
g/L  cellulose  was  used.   When  instability  of  the  reactor  was 
observed,  the  feedstock  concentration  was  reduced  to  22.5 
g/L  and  kept  constant  for  the  rest  of  the  experiment. 

During  steady  state  operation  in  the  first  phase,  three 
different  flow  rates  (or  3  different  organic  loading  rates, 
since  the  feedstock  concentration  was  constant)  were  used. 
Performance  data  collected  in  this  stage  were  analyzed  and 
used  to  initialize  the  dynamic  control  operation  for  the 
second  phase. 
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Three  tests  were  conducted  for  the  second  phase.   One 
was  for  the  objective  of  "maximum"  methane  production  rate 
(MPR)  operation  and  another  two  were  designed  for  the 
designated  MPR  operations  at  set  points  of  0.4  L  CH4/L-D  and 
0.2  L  CH4/L-D.   The  set  point  of  0.4  L  CH4/L-D  was  selected 
because  it  was  about  80%  of  the  maximum  MPR  achieved  in  the 
steady  state  operation.   The  set  point  of  0.2  L  CH4/L-D  was 
chosen  to  test  the  flexibility  of  the  system,  i.e.,  to 
examine  whether  it  would  function  at  both  high  (0.4)  and  low 
(0.2)  set  points.   The  "maximum"  MPR  was  defined  as  the 
maximum  MPR  which  could  be  achieved  under  the  current 
reactor  conditions  without  causing  any  instability  of  the 
reactor.   The  operating  procedures  for  each  test  are 
summarized  in  Table  3-4. 


Table  3-4.   Operating  procedures  for  the  anaerobic 
digestion  using  cellulose  as  feedstock. 


Test  no.  Operation 


First  phase:  (steady  state) 


a 


51  at  a  flow  rate  of  250  ml/day 

52  at  a  flow  rate  of  375  ml/day 

53  at  a  flow  rate  of  450  ml/day 

Second  phase:  (dynamic  control) 

D*  set  point  MPRsp  =  0.4  L  CH4/L-day 

D2         maximum  MPR 

D3         set  point  MPRsp  =  0.2  L  CH4/L-day 

a.  "S"  indicates  the  steady  state  operation. 

b.  "D"  indicates  the  dynamic  control  operation. 
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The  reactor  was  fed  every  12  hours  (twice  a  day)  for 
both  phases.   The  sampling  frequency  was  every  two  or  three 
days  in  the  first  phase  and  twice  a  day  for  the  second 
phase.   Chemical  oxygen  demand  (COD) ,  total  organic  carbon 
(TOC) ,  pH,  gas  production  and  gas  composition  were  measured 
for  each  sample.   Total  solids  (TS) ,  volatile  solids  (VS) 
and  volatile  acids  (VA)  were  monitored  the  same  frequency  as 
other  measurements  in  the  first  phase  and  were  changed  to 
once  every  three  days  for  the  second  phase.   Because  the 
determination  of  COD  takes  much  longer  than  that  of  TOC  (2-3 
hours  compared  with  15-20  minutes) ,  COD  was  estimated  from 
correlations  with  TOC  during  dynamic  control  in  the  second 
phase.   The  sampling  frequencies  and  variables  monitored  for 
both  phases  are  summarized  in  Table  3-5. 

Analytical  Methods 

Total  solids  (TS)  and  volatile  solids  (VS)  were 
analyzed  according  to  the  procedures  of  Standard  Methods 
(APHA,  1985) .   Total  volatile  acids  concentrations  were 
analyzed  using  the  methods  of  DiLallo  and  Albertson  (1961). 
pH  was  measured  with  an  Orion  701  digital  ionalyzer. 

Chemical  oxygen  demand  (COD)  was  determined  using  Hach 
COD  vials  (Hach  Co.,  Colorado).   The  titration  method  was 
used  at  the  beginning  of  the  experiment,  and  the 
colorimetric  method  was  used  after  comparing  two  lots  of 
standard  solution  for  both  methods.   The  results  showed  that 
the  colorimetric  method  was  more  accurate  than  the  titration 
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Table  3-5.   Sampling  frequency  and  measured  variables  for 
the  anaerobic  digestion  using  cellulose  as 
feedstock. 


Phase  1 


Phase  2 


Sampling  frequency 
Measured  variables3 

TS,    g/L 

VS,    g/L 

COD,    g/L 

TOC,    g/L 

VA,    g/L 

PH 

Gas  production,  L 

Gas  composition 


every  2  or  3  days   every  12  hours 


Yes 
Yes 
Yes 
Yes 
Yes 
Yes 
Yes 
Yes 


YesL 

Yesfc 

Yes 

Yes 

Yesb 

Yes 

Yes 

Yes 


a.  Measurements  were  conducted  for  each  sample  except 
those  indicated. 

b.  Measured  once  every  3  days. 
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method.   Transmittance  against  distilled  water  was  read  at 
620  nut  on  a  Bausch  &  Lomb  spectrophotometer  (Spectronic  70)  . 
Standard  curve  and  regression  line  were  calculated  by  using 
the  standard  solution  prepared  from  potassium  hydrogen 
phthalate  (KHP)  for  each  set  of  samples. 

Total  organic  carbon  (TOC)  was  measured  by  using  an 
Astro  1200  Carbon  Analyzer  combined  with  an  Astro  5000 
Infrared  Analyzer  (Astro  Resources  Corp.,  Texas).   The 
sample  was  acidified  by  using  a  solution  of  5  N  Phosphoric 
acid  (H3P04)  to  bring  the  sample  pH  within  a  range  of  2  to  3 
before  injecting  into  the  analyzer.   Standard  curve  (Figure 
3-3)  of  TOC  concentration  versus  analog  output  (millivolt) 
was  made  by  using  the  standard  solution  prepared  from 
potassium  hydrogen  phthalate  (as  stated  in  Standard  Methods, 
APHA,  1985) .   To  assure  its  reliability,  the  standard  curve 
was  checked  and  modified  (if  necessary)  each  week.   The 
first  5  minutes  after  injection  for  each  new  sample  was 
skipped  as  the  stabilization  time  and  the  analog  output  was 
read  every  10  seconds  for  a  period  of  10  minutes.   A  Zenith 
248  microcomputer  and  a  Keithley  570  workstation  were 
combined  to  receive  analog  outputs  (as  mV)  from  the  infrared 
analyzer  and  display  the  TOC  concentration  (as  mg/L)  on  the 
console. 

Gas  production  was  measured  before  each  sampling  with  a 
Wet  Test  Gas  Meter  (Model  63115,  Precision  Scientific, 
Chicago,  Illinois)  at  a  room  temperature  of  27  °C.   Gas 
composition  was  determined  on  a  Gow-Mac  550  Gas 
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Chromatograph  with  a  thermal  conductivity  detector  (Gow-Mac 
Inst.  Co.,  Madison,  New  Jersey).   A  stainless  steel  column 
(2.44  m  by  6.35  mm)  packed  with  50/80  Porapak  Q  (Supelco, 
Bellefonte,  Pennsylvania)  was  used.   Helium  was  used  as  the 
carrier  gas  at  a  flow  of  60  ml/min.   The  temperatures  of 
injection  port,  column  and  detector  were  90,  70  and  150  °C, 
respectively.   Results  were  recorded  and  integrated  by  using 
an  Apple  II  compatible  computer  which  was  eguipped  with  a 
Chromcard  (Anadata,  Inc.,  Glen  Ellyn,  Illinois)  for 
chromatography  data  handling. 

Controller 

The  controller  is  an  element  which  receives  the 
information  from  the  measuring  devices,  makes  decisions 
based  upon  the  control  algorithm,  and  takes  the  appropriate 
control  action  to  adjust  the  values  of  the  manipulated 
variables.   Therefore,  the  controller  in  a  computer  control 
system  is  a  combination  of  hardware  and  software.   This 
section  will  discuss  only  the  software  related  to  the 
control  algorithm,  and  the  hardware  and  the  control  action 
(or  actuating  signal)  generating  program  will  be  described 
in  a  later  section. 

Model  Development 

Because  it  has  well  defined  characteristics,  a 
continuous  stirred  tank  reactor  (CSTR)  was  selected  for  this 
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control  system.   The  mathematical  model  is  a  set  of  state 
equations  to  relate  the  state  variables  to  the  various 
independent  variables.   The  principle  of  conservation  states 
that: 

Accumulation  =  Input  -  Output  +  formation  -  consumption 

The  mathematical  model  of  this  anaerobic  digestion  system 
was  developed  based  on  the  mass  balance  of  the  system 
components  shown  in  Figure  3-4.    Although  four  groups  of 
bacteria  have  been  recognized  in  anaerobic  digestion 
processes,  two-culture  Monod  growth  kinetics  were  adopted 
here  because  of  their  simplicity  and  accuracy  for  real-time 
control.   The  bioconversion  steps  from  cellulose  to  biogas 
are  shown  below: 


Cellulose  — - — >  Volatile  acids  +  CO,  +  Cells 

Acetogens  z 


Volatile  acids  „  ., >  CO,  +  CH,  +  Cells 

Methanogens     2     4 


Definition  of  variables  and  parameters 

The  definition  of  variables  (including  input,  output 
and  state  variables)  and  parameters  used  in  this  study  are 
listed  in  Tables  3-6  and  3-7,  respectively.   A  complete 
listing  of  notations  is  shown  in  Appendix  A. 


38 


CH4+  CO  2 


Influent 


Q.Si.Vai.Xai.Xmi 


Q,S,Va,Xa,Xm 


Figure  3-4.   Schematic  diagram  of  the  CSTR  anaerobic 
digestion  system. 
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Table  3-6.  Definition  of  variables  used  in  the  control 
system. 


Variable 


Definition 


s 

Va 
Xa 
Xm 

so 
Va0 

Xa0 

Xm0 

Si 

Vai 

Xai 

Xmi 

Q 

V 

**1 

M2 

QCH4 
QC02 

Qgas 
GPR 

MPR 

MPR 

MPR 


sp 


Substrate  concentration  in  reactor,  g/L 

Total  volatile  acids  cone,  in  reactor,  g/L 

Acetogen  concentration  in  reactor,  g/L 

Methanogen  concentration  in  reactor,  g/L 

Initial  substrate  concentration  in  reactor,  g/L 

Initial  total  volatile  acids  cone,  in  reactor,  g/L 

Initial  acetogen  concentration  in  reactor,  g/L 

Initial  methanogen  concentration  in  reactor,  g/L 

Influent  substrate  concentration,  g/L 

Influent  total  volatile  acids  cone. ,  g/L 

Influent  acetogen  concentration,  g/L 

Influent  methanogen  concentration,  g/L 

Flow  rate,  L/day 

Manipulated  flow  rate,  L/day 

Liguid  volume  of  reactor,  L 

Specific  growth  rate  of  Xa,  l/day 

Specific  growth  rate  of  Xm,  l/day 

Daily  methane  production,  L  CH4/day 

Daily  carbon  dioxide  production,  L  C02/day 

Daily  biogas  production,  L/day 

Volumetric  biogas  production  rate,  L/L-day 

Volumetric  methane  production  rate,  L  CH4/L-day 

Predicted  methane  production  rate,  L  CH,/L-day 

Set  point  of  volumetric  methane  production  rate, 

L  CH4/L-day 


Table  3-7.  Definition  of  parameters  used  in  the  control 
system. 


Parameter 


Definition 


MB1 

Mm2 
Ks1 

Ks2 

Kd1 

Kd2 

Ya 

Yb 

Yc 

Yd 

Ye 

Yf 


Maximum  specific 
Maximum  specific 
Saturation  coeff 
Saturation  coeff 
Decay  rate  coeff 
Decay  rate  coeff 
Xa  yield  coeff. , 
Va  yield  coeff. , 
Xm  yield  coeff. , 
CH4  yield  coeff. , 
C02  yield  coeff. 
C02  yield  coeff. 


growth  rate  of  Xa,  day"1 

growth  rate  of  Xm,  day"1 
.  for  Xa,  g/L 
.  for  Xm,  g/L 
.  for  Xa,  day"1 
.  for  Xm,  day"1 

g  Xa  prod./g  S  used 

g  Va  prod./g  S  used 

g  Xm  prod/g  Va  used 
L  CH4  prod./g  Va  used 
via  acetogen,  L  C02/g  Xa  prod, 
via  methanogen,  L  C02/g  Xm  prod. 
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Mass  balance 

Substrate  (S) .   Substrate  used  for  the  mathematical 
model  can  be  any  indicator,  such  as  volatile  solids, 
chemical  oxygen  demand  or  organic  carbon,  which  is  able  to 
reflect  the  organic  level  of  the  digester. 


Change  of  S  =  S  from  the  influent  -  S  in  the  effluent 
-  utilization  of  S  to  form  acetogen  (Xa) 


dS  jLt^Xa-V 

V-( )  =  Q-Si  -  Q«S (3-1) 

dt  Ya 


where  /j1  = 


Ks^S 


Volatile  acids  (Va) . 

Change  of  Va  =  Va  from  the  influent  -  Va  in  the 
effluent  +  formation  from  S  - 
utilization  for  methanogen  (Xm)  growth 


dVa  Yb-^'Xa'V     M2*Xm*v 

V»  ( )  =  Q-Vai  -  Q«Va  + 


dt  Ya  Yc 

(3-2) 

where  /x2  =  

Ks2+S 


Acetogen  (Xa) . 

Change  of  Xa  =  Xa  in  influent  -  Xa  in  effluent  +  Xa 

formed  from  S  -  decay  of  Xa 


dXa 

V*  ( )  =  Q«Xai  -  Q«Xa  +  /^'Xa-V  -  Kd^Xa'V     (3-3) 

dt 
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Methanoqen  (Xm) . 

Change  of  Xm  =  Xm  in  influent  -  Xm  in  effluent  +  Xm 

formed  from  Va  -  decay  of  Xm 


dXm 

V- ( )  =  Q'Xmi  -  Q-Xrn  +  n2'Xm-V   -  Kd2'Xm»V    (3-4) 

dt 


Gas  production.   Methane  and  carbon  dioxide  production 
were  calculated  based  on  the  cell  mass  produced,  and  the  gas 
production  was  the  sum  of  methane  and  carbon  dioxide 
produced. 


V'Yd-Mz'Xm 

>CH4  =  (3-5) 

Yc 


Qco2  =  Ye-^'Xa  +  Yf»/i2'Xm  (3-6) 

Qgas   "   Qch4   +    Qco2  (3-7) 

Volumetric  methane  production  rate  (MPR)  was  found  by 
integrating  the  methane  produced  during  a  certain  period  and 
then  divided  by  the  time  length  of  that  period. 


j; 


t+<St 

QCH4dt 


MPR   =   (3-8) 

V«5t 

Volumetric  gas  production  rate  (GPR)  was  found  in  the  same 

way  as  methane  production  rate  but  including  both  methane 

and  carbon  dioxide  produced. 

>t+(5t 

QGASdt 

GPR   = (3_9) 

V-6t 
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The  above  model  can  be  used  to  simulate  the  dynamic 
behavior  of  the  digester  and  find  the  desired  variables  with 
the  entry  of  parameter  values  and  initial  conditions. 
Identification  of  the  required  parameters  and  the  model 
validation  will  be  discussed  in  a  later  section. 

Control  Algorithm 

Selecting  the  proper  measurements  and  the  variable  for 
control  is  essential  for  establishing  a  control  algorithm. 
Gas  production,  gas  composition,  and  total  organic  carbon 
(TOC)  were  selected  as  the  measured  variables  for  the 
control  system.   Gas  production  and  gas  composition  were 
combined  to  find  the  methane  production  rate — the  output  (or 
controlled)  variable.   Total  organic  carbon  was  chosen 
because  of  its  on-line  characteristics  and  being  directly 
related  to  the  organic  matter  content.   The  flow  rate  (Q) 
was  selected  as  the  sole  manipulated  variable  because  of  its 
ease  of  operation  and  close  relationship  with  the  organic 
loading  and  the  methane  production  for  a  given  feedstock. 

A  variety  of  algorithms,  from  the  conventional 
proportional  control  to  the  modern  adaptive  control,  are 
available  for  process  control  (Stephanopoulos,  1984;  Banks, 
1986;  Balchen  and  Mumme,  1988;  Landau,  1979;  Chalam,  1987). 
To  select  an  appropriate  control  algorithm  for  the  system, 
both  the  mathematical  model  and  the  objective  function  need 
be  considered.   The  objectives  for  this  study  were  to 
operate  the  anaerobic  reactor  (1)  at  a  designated  methane 
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production  rate  (minimize  |MPRsp-MPR|)  and,  (2)  at  the 
maximum  methane  production  rate  (maximize  MPR)  without 
causing  the  failure  of  the  reactor.   In  order  to  maintain 
reactor  stability,  a  constraint  of  minimum  organic  removal 
was  added  to  the  objective  function.   Therefore,  the  control 
algorithm  for  this  system  had  to  conduct  a  set-point  control 
as  well  as  to  predict  the  state  variables  for  the  system. 
The  model  proposed  above  is  a  nonlinear  model.   The 
preliminary  results  showed  that  this  system  was  unobservable 
which  meant  that  the  state  variables  could  not  be  completely 
observed  from  the  measurements  of  the  outputs  (Kuo,  1987; 
Banks,  1986) .   Since  the  prediction  of  the  state  variables 
and  the  set-point  variation  were  requisite  for  this  control 
system,  a  nonlinear  self -tuning  regulator  was  chosen  to 
achieve  the  objective  function. 

Nonlinear  self-tuning  regulator  (NSTR) 

The  nonlinear  self-tuning  regulator  (NSTR) ,  a  type  of 
adaptive  controller  (Stephanopoulos,  1984;  Chalam,  1987),  is 
the  combination  of  an  adaptive  parameter  estimator  (PE)  and 
an  optimizer  (OP).   As  shown  in  Figures  3-5  and  3-6,  the 
database  required  for  running  PE  and  OP  is  updated  with  the 
measured  outputs  (including  the  controlled  output  MPR)  and 
the  calculated  control  variable  (Qm,  for  current  feeding)  at 
each  sampling  time.   Using  the  least-square  principle,  a 
parameter  estimator  generates  a  set  of  parameters  to  fit  the 
historical  data.   The  optimizer  then  revises  the  simulation 
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model  with  the  above  parameters  and  finds  the  optimum 
manipulated  variable  (Qm,  for  next  feeding)  for  the  next 
sampling  time  to  approach  the  desired  output  (MPR   or 
maximum  MPR  )  . 

Parameter  estimator 

Method.   The  algorithm  used  in  this  study  for  parameter 
estimation  is  modified  from  the  "Complex"  method  developed 
by  Box  (1965) .   The  Complex  method  is  a  sequential  search 
technique  for  solving  problems  with  a  nonlinear  function  of 
multiple  variables  subject  to  nonlinear  inequality 
constraints.  In  other  words,  it  searches  for  the  objective 
function  in  a  feasible  region  defined  by  the  upper  and  lower 
bounds.   This  method  requires  no  derivatives  and  should  tend 
to  find  the  global  optimum  because  the  initial  set  of  points 
are  randomly  scattered  throughout  the  feasible  region.   In 
the  modified  method  used  in  this  study,  the  objective 
function  was  set  to  minimize  total  residual  sum  of  squares 
of  GPR  (gas  production  rate) ,  MPR  (methane  production  rate) , 
PCH4  (methane  content)  and  COD  (chemical  oxygen  demand) . 
The  values  of  parameters  were  adjusted  recursively  inside 
the  defined  region  to  approach  the  objective  function.   The 
residue  was  found  by  comparing  the  simulation  results  and 
the  observed  data.   Therefore,  the  set  of  parameters  which 
satisfied  the  objective  function  fitted  the  curve  the  best. 
To  avoid  the  bias  caused  by  the  different  scale  between  each 
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variable,  the  residue  was  corrected  as  the  fraction  of  the 
observed  data  as  follows: 


Y  -  Yp 
corrected  residue  =  


where  Y  is  the  observed  data  and  Yp  is  the  predicted  data 
from  simulation. 

Range  for  searching  parameters.   In  order  that  the 
parameters  found  have  physical  meaning,  the  upper  and  lower 
bounds  were  carefully  selected  by  referring  to  several 
kinetic  studies  of  the  anaerobic  digestion  process.   The 
parameter  values  from  the  referred  studies  are  listed  in 
Table  3-8,  and  the  bounds  of  search  were  initially  set  by 
using  or  covering  the  extreme  points  from  these  studies. 
However,  preliminary  results  showed  that  error  code  occurred 
for  some  parameter  values  during  execution  of  the  program. 
Therefore,  modifications  were  made  for  those  parameters  to 
successfully  execute  the  computer  program.   The  final  ranges 
for  searching  are  listed  in  Table  3-9.   The  initial 
conditions  for  Xa  and  Xm  (Xa0  and  Xm0) ,  which  were  unknown, 
were  added  to  the  estimation  process  even  though  they  were 
not  kinetic  parameters.   The  searching  bounds  for  Xa0  and 
Xm0  were  set  as  1%  -  20%  and  0.1%  -  10%,  respectively,  of 
the  initial  VS  concentration  of  the  reactor. 
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Table  3-9.   Range  of  search  for  parameter  values. 


Parameter        Lower  bound    Upper  bound 

Ya,  g  Xa/g  S  0.1  1.0 

Yb,  g  Va/g  S  0.05  2.0 

Yc,  g  Xm/g  Va  0.03  1.0 

Yd,  L  CH4/g  Va  0.3  5.0 

Ye,  L  C02/g  Xa  0.5  2.0 

Yf,  L  C02/g  Xm  2.0  10.0 

Ks1(  g/L  0.1  10.0 

Ks2,  g/L  0.15  10.0 

Kd1#  day"1  0.01  0.5 

Kd2,  day"1  0.01  0.05 

/im1,  day"1  0.04 


5.0 


Mncr    day"4  0.004  1.5 

Xa0,    g/L#  0.1 


2.0 


Xm0,  g/L*  0.01  1.0 

*  Variable  included  in  the  estimation  procedure 
but  not  covered  in  Table  3-8. 


Algorithm.   The  parameter  searching  algorithm  proceeded 

as  follows: 

1.  A  feasible  starting  point  G1  ,  (a  set  of  model  parameters 
of  which  the  parameter  values  were  within  the  range 
listed  in  Table  3-9)  was  chosen,  and  N  additional  points 
were  generated  by  using  the  pseudorandom  numbers  r,  ,  and 

*  r  J 

the  constraints  for  each  individual  parameters  (MAXIM- 
and  MINIMj)  to  make  an  N+l  by  N  matrix  as  the  initial 
guesses  for  the  parameters. 

Gi.j  =  {G1,J  +  [MINIMj+r,-^.' (MAXIMj -MINIMj)  ]}  +  2     (3-10) 
i  =  2, N+l  and  j  =  i,N 

where  Gi(j  is  the  jth  parameter  (same  order  as  in  Table  3- 
9)  of  the  ith  point  (each  point  is  a  set  of  parameters) , 
N  is  the  number  of  parameters  to  be  estimated  (N  is  equal 
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to  14  in  this  study) ,  G. ,  is  the  feasible  starting  point, 
MINIM-  and  MAXIM-  are  values  of  the  jth  lower  and  upper 
bound  as  in  Table  3-9,  respectively,  and  r,  ,  are 

'  i  J 

uniformly  distributed  over  the  interval  [0,1].   This 
algorithm  requires  N+l  points  to  estimate  N  parameters. 

2.  The  objective  function,  total  residual  sum  of  squares 
(RSS) ,  was  evaluated  at  each  point  based  on  the  initial 
guesses  found  above: 

RSS  =  RSSGPR  +  RSSMPR  +  RSSCOD  +  RSSPCH4 
where  RSSGPR,  RSSMPR,  RSSCOD  and  RSSPCH4  are  the  residual 
sum  of  squares  for  GPR,  MPR,  COD  and  methane  content 
(PCH4) ,  respectively.   Residual  sum  of  squares  for  each 
variable  were  calculated  as  the  sum  of  the  corrected 
residues  (as  described  in  the  previous  section)  during 
the  simulation  period.   Within  the  N+l  RSS's,  the  best 
guess  (minimum  RSS)  and  worst  guess  (maximum  RSS)  were 
selected  and  assigned  indexes  as  "ic"  and  "ir" , 
respectively . 

3.  A  new  set  of  parameters  (new  point)  was  determined  and 
added  to  the  bottom  of  the  parameter  matrix  as  follows: 

N+1 

(l+a)..S   Gj  , 

1=1  J  ' 

GN+2.i    =  n f-N        +    1>'6lp.i  (3"11> 

a   >   1   and   i  =   1,N 

where  Gjr  {  was  the  set  of  parameters  which  produces  the 
maximum  RSS  and  a  was  the  over-reflection  factor.   Box 
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(1965)  recommended  a  value  of  a=1.3,  and  a=l  was  selected 
and  used  here  after  some  preliminary  runs. 

4 .  RSS  of  the  new  guess  was  evaluated  and  compared  with  RSS 
of  the  worst  guess.   The  worst  guess  was  replaced  if  its 
RSS  was  greater  than  that  of  the  new  guess  (Gir  ,-  =  GN+2  ,-, 
i=l,N) .   If  RSS  of  the  new  guess  was  greater  than  or 
egual  to  that  of  the  worst  guess,  each  point  was  moved 
one  half  of  the  distance  to  the  best  guess. 

Gi.j  =  <Gi.J  +  Gic,j)  +  2  (3-12) 

where  Gjc  .  was  the  set  of  parameters  which  produced  the 
minimum  RSS. 

5.  All  the  selected  points  had  to  be  located  in  the  feasible 
region.   In  case  that  a  parameter  for  a  guess  was  beyond 
this  region,  the  value  of  the  lower  or  upper  bound  was 
assigned  to  that  particular  parameter,  depending  on  which 
bound  was  closer. 

6.  A  stopping  rule  was  set  at  the  time  at  which  the 
parameter  matrix  was  reduced  to  an  acceptably  small  size 
e,  (convergence  test)  or  the  RSS  of  the  best  guess  was 
less  than  a  predetermined  small  number  e2  (minimum  RSS 
test) .   The  convergence  test  was  conducted  as  follows: 

{^h~l!  tRSSc  "  RSS,.]2}1/2  <  e,  (3-13) 

where  RSSC  is  the  RSS  of  the  centroid  and  RSS,.  is  the  RSS 
of  the  ith  guess.   The  centroid  is  determined  by: 
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CENTR0j  =  -£-  (  "£   Gjl.  -  Gfrf|)     ,i  =  1,N      (3-14) 


where  CENTRCK  was  the  ith  parameter  of  the  centroid. 

The  flow  chart  of  the  parameter  estimator  is  shown  in 
Figure  3-7  and  the  program  is  listed  in  Appendix  B. 

Optimizer 

Method.   Fibonacci  search  method  was  used  to  find  the 
optimum  value  of  the  control  variable  (flow  rate,  Qm)  .   It 
is  an  interval  elimination  search  method  for  the  unimodal 
function  (Jacoby  et  al.,  1972;  Dixon,  1972;  Kuester  and 
Mize,  1973) .   The  advantages  of  the  Fibonacci  method  are  no 
requirement  of  derivatives,  guarantee  of  a  prescribed 
accuracy,  and  simplicity  of  coding. 

Similar  to  the  parameter  estimator,  the  optimizer  also 
used  the  simulation  results  to  calculate  the  values  of 
objective  function.   However,  the  parameter  estimator  only 
simulated  the  real  experimental  data,  whereas  the  optimizer 
extended  its  simulation  one  more  step  to  the  next  sampling 
time  which  had  no  real  data  with  which  to  compare. 
Depending  upon  the  purpose,  the  objective  functions  were  set 
to  minimize  the  square  residue 

Min.  (MPRsp  -  MPRp)2, 

or  to  maximize  the  controlled  output 


Max.  MPRp. 


(    START    ) 
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INPUT  FEASIBLE  STARTING  POINT,  G  i,i    ,j=1,N 


GENERATE  ANOTHER  N  SETS  OF  PARAMETERS: 
Gjj  =   (G)ij+   (MINIMj  +  qj  (MAXIM  j- MINIM  j)))/2  ,  i=2,N,  j=1,N 


■> 


CALL  SIMULATION   MODEL  AND  FIND  RSS  FOR  EACH  POINT 


DETERMINE  THE  BEST  (IC)  AND  THE  WORST  (IR)  GUESS 


YES 


YES 


DETERMINE  "NEW"   SET  OF  PARAMETERS 


(1  +  0)51  G  1 1 

GN+2.i  = FT-^-^  +  QGw  ■    a*',  i=i,N 


YES 


">1  GN+2,i    =  MAXIMj 


YES 


~H       GN-f2,i  =  MINIMj" 


CALL  SIMULATION  MODEL  AND  FIND  RSS(N  +  2)  FOR  THE  "NEW  GUESS 


MOVE  ONE  HALF  OF 
THE  DISTANCE  TO 
THE  BEST  GUESS 
Gi,j=  (Gi.j+GiC.j)/2 


(    STOP     ) 


Figure  3-7.   Flow  chart  of  the  parameter  estimator. 
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Both  functions  were  subject  to  the  following  constraints: 


CODn  <  CODc„  and 

P        sp 


0  <  Q  <  QHIGH 


-m 


where  COD  is  the  predicted  COD  from  the  simulation  model, 
COD   is  the  predetermined  COD  concentration,  Qm  is  the 
manipulated  flow  rate  and  QHIGH  is  the  upper  bound  of  the 
flow  rate. 

The  COD  constraint  was  set  to  maintain  a  minimum  limit 
of  85%  organic  removal  efficiency  to  avoid  failure  of  the 
process.   Selection  of  the  original  searching  region  was 
critical  for  the  final  solution,  and  it  was  possible  that 
there  were  more  than  one  value  that  could  satisfy  the 
objective  function.   However,  a  conservative  range  of  flow 
rate  from  0  to  1  L/day  (i.e.  HRT  >  10  days)  was  selected  and 
runs  with  different  maximum  flows  (0.25,  0.5,  0.75  and  1.0 
L/day)  were  made  to  circumvent  these  problems.   In  the  case 
of  multiple  solutions,  the  smallest  value  was  used  for  safe 
operation. 

Algorithm.   The  algorithm  for  designated  MPR  operation 
proceeded  as  follows: 

1.  An  original  search  region  R,  was  specified  with  QLOW  and 
QHIGH  as  the  lower  and  upper  boundaries,  respectively, 
where 

R,    =  QHIGH  -  QLOW.  (3-15) 

2.  The  predetermined  desired  accuracy  a  was  defined  as  the 
ratio  of  the  final  region  divided  by  the  original  search 
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region.   Using  the  value  of  a,    the  largest  Fibonacci 
number  FN  was  defined  as 

a  <  — I (3-16) 

and  F0  =  F1  =  1, 

Fn  =  Fn  ,  +  Fn  ,,  for  n  >  2 

n      n-l      r\-d' 

where  Fn  is  called  a  Fibonacci  number,  and  N  is  egual  to 
the  number  of  iterations  of  the  search  procedure. 
3.  Find  the  first  two  interior  points  Q1  and  Q2  (Q.,  <  Q2)  , 
where 


p 

Q,  =  QLOW  + R,  (3-17) 

F 


F 

Q2  =  QHIGH — R,  (3-18) 

F 

and  compute  values  of  the  objective  function,  M(Q^)    and 
M(Q2)  ,  at  these  two  points.   M^)  and  M(Q2)  were 
calculated  as  the  sguare  residues  for  MPR  at  different 
flow  rates  Q,  and  Q2,  respectively. 

|MPRp  -  MPRsp|2 
where  MPR  is  the  predicted  MPR  and  MPRCP,  is  the  desired 
MPR. 
4.  Adjust  the  search  region  as  follows: 

If  M(Q2)  <  M(Q,)  ,  the  lower  region  [QLOW^)  was 
discarded,  and  the  boundary  was  reset  as  QLOW  =  Q, . 
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The  new  search  region  was  defined  as 


*Vi 

FN-2 

•  R1 

i<2    —                                   •     K^     —    K^     — 

FN 

F 

while  the   interior  points, 

Q1    =   Q2,    and 

FN-3 

^2  ~~  yiixoii   —                          •    Kp . 

FN-1 

(3-19) 


(3-20) 


If  M(Q2)  >  M(Q1)  ,  the  upper  region  (Q2,QHIGH]  was 
discarded,  and  the  boundary  was  reset  as  QHIGH  =  Q2. 
The  new  search  region  was  the  same  as  the  expression  in 
Eg.  3-19,  and  the  interior  points  were  defined  as 
follows: 


R2,  (3-21) 


FN-3 

y1   —  yijuw   + 

FN-1 

Q2   =   Q,. 

5.  Continue  the  above  region  reducing  procedure  for  N 
evaluations  using  the  following  general  eguations: 


FN-(k-1)  FN-k 


Rk  = R,  =  Rk_, R,.,       (3-22) 

FN 

and  the  interior  points 


FN  FN-(k-2) 


FN-(k+1) 
Q,  =  QLOW  +  •  Rk  (3-23) 

FN-(k-1) 
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or 

FN-(k+1) 
Q2  =  QHIGH •  Rk.  (3-24) 

FN-(k-1) 

6.  After  discarding  the  proper  region  in  each  iteration,  the 
last  two  test  points  would  be  coincident  and  located  in 
the  center  of  the  remaining  region.   To  determine  which 
half  of  the  remaining  region  the  optimum  belonged  to,  one 
of  the  last  two  points  was  shifted  by  a  small  distance  e. 
The  objective  function  was  then  evaluated  at  this  point 
and  the  final  region  where  the  optimum  was  located  would 
be  determined.   The  location  of  the  optimum  point  was 
assumed  to  be  the  midpoint  of  the  final  region. 

7.  Check  if  the  COD  concentration  for  the  optimum  point  was 
within  the  limit  of  the  constraint.   If  the  COD  was 
beyond  the  specified  limit,  the  value  of  the  optimum 
point  was  reduced  by  0.05  L/day  each  time  until  the 
resulting  COD  was  within  the  limit.   The  objective 
function  and  the  predicted  MPR  were  also  evaluated  for 
this  "final"  optimum  point. 

Besides  the  difference  in  the  objective  function,  the 
algorithmic  structure  of  the  maximum  MPR  operation  was  the 
same  as  the  above  designated  MPR  operation.   The  only 
modification  reguired  was  inverting  the  "<"  to  ">"  in  the 
"IF"  statement  of  step  (4) . 

The  flow  chart  of  the  optimizer  is  shown  in  Figure  3-8 
and  the  program  listing  for  the  optimizer  is  shown  in 
Appendix  B. 


(     START    ) 
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INPUT  VALUES  OF  PARAMETERS  FOUND  IN  PE 


DEFINE  THE  DESIRED  ACCURACY  «  AND  THE 
ORIGINAL  SEARCH  REGION  R,,  QHIGH,  QLOW 


Ri+1=  QHIGH   -   QLOW 

Q2  =  Q1 
i 


Q1    =  QLOW  +   F|k     *  Ri+i/Fjj  Q2   =  QHIGH   -   F,k     *  Ri+1/Fjj 


Ri+i=  QHIGH  -  QLOW 

,  t  i 

Q1    =  Q2 


X 


EVALUATE  OBJECTIVE  FUNCTION   M(Q1)  AND  M(Q2) 


ik=ik— 1,  jj=jj-1,  i=i+1 


£     =  0.001   *  Q1 


DL  =  Q1   +    £ 


EVALUATE  OBJECTIVE  FUNCTION  M(DL) 


NO 


QM  =  (Q1    +QL0W)/2 

QM  =  (Q2  +  QHIGHV2 

I 

■ 

f 

ACCURACY  =  (Q1-QL0W)/R, 

>- 

LZ 

EVALUATE  OBJECTIVE  FUNCTION  M(QM)     I 

NO 

=CTIC0D^C0Di  )/COD;  ^J55?Tl_t>- 

1 

QM  =  QM  -  0.05 

^"J'yes 

I  OUTPUT  QM| 

(   STOP    ) 


Figure  3-8.   Flow  chart  of  the  optimizer. 
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Hardware 

In  order  to  verify  the  above  proposed  control 
algorithm,  a  real-time  control  system  was  set  up  and 
different  operational  objectives  were  tested.   A  schematic 
diagram  of  the  integrated  control  system  is  shown  in  Figure 
3-9.   For  the  designated  MPR  operation,  0.2  and  0.4  L  CH4/L- 
day  were  selected  as  the  set  points.   The  maximum  MPR 
operation  was  also  tested. 

This  section  will  cover  the  hardware  used  in  the 
control  system  and  the  required  programming  to  generate  the 
control  action.   The  hardware  consisted  of  computer,  data 
acquisition  system,  analog-digital  (A/D  and  D/A)  converters 
and  the  final  control  elements. 

Computer 

Computers  were  used  in  this  system  for  input  of  the  TOC 
(via  an  A/D  converter) ,  estimating  model  parameters,  finding 
the  optimum  value  for  the  control  variable  and  sending 
digital  signals  (turn  flow  pumps  on  or  off)  to  take  the 
control  action.   In  order  not  to  interrupt  each  other,  two 
different  computers  were  used.   As  shown  in  Figure  3-9,  both 
parameter  estimator  and  optimizer  were  coded  in  FORTRAN  and 
run  in  a  Sun  3/260  minicomputer  where  the  database  was 
stored.   The  reason  for  choosing  the  Sun  computer  was  its 
rapid  speed  of  calculation.   The  other  computer  used  in  this 
system  was  a  Zenith  248  microcomputer  which  served  as  a 
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61 
supervisor  to  receive  the  information  from  the  optimizer  and 
then  generate  the  proper  control  action  at  the  right  time. 
To  accomplish  the  supervisory  function,  a  control  program 
was  developed  to  calculate  the  required  pumping  time  based 
on  the  optimum  flow  rate  found  in  the  optimizer.   The 
control  program  also  performed  as  a  timer  and  sent  the 
signals  to  the  following  power  control  system.   To  work  as  a 
timer,  this  control  program  had  a  loop  to  read  the  computer 
clock  for  continuously  monitoring  the  process.   This  control 
program  was  coded  in  a  programing  language,  SOFT500,  which 
is  unigue  to  the  Keithley  570  system  and  compatible  with  the 
BASIC  programming  language.    The  flow  chart  for  the  control 
program  is  shown  in  Figure  3-10,  and  the  program  code  is 
listed  in  Appendix  B. 

Data  Acquisition  Workstation 

The  data  acquisition  system  used  to  receive  signals 
from  the  microcomputer  was  a  Keithley  570  workstation 
(Keithley  Data  Acquisition  &  Control,  Cleveland,  Ohio)  which 
contained  different  slots  to  achieve  functions  of  analog 
I/O,  digital  I/O  and  relay  control.   The  communication 
between  the  microcomputer  and  the  workstation  was 
accomplished  by  an  interface  card  which  was  plugged  into  a 
full  length  slot  of  the  microcomputer.   The  signal  sent  to 
the  workstation  was  then  transmitted  to  a  power  control 
board  (Opto  22,  Huntington  Beach,  California)  which 
contained  plug-in  relays  to  regulate  the  on/off  of  the  flow 
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("    START  ~) 


Y 

INPUT  THE  OPTIMUM   FLOW  RATE,   QM 


"■ 


CALCULATE  THE  REQUIRED  PUMPING  TIME 


RESET  THE  BINARY  CODE  OF  THE  RELAY: 
SWITCH   =  0 


DEFINE  A  VARIABLE  "RELAY"   ON  THE 
RELAY  CONTROL  SLOT 


NO 


YES 


■■-' 


CHANGE  THE   BINARY  CODE  TO   "ON"  VALUE 
SWITCH   =   2 


"' 


ASSIGN  THE  NEW  VALUE  TO  THE  VARIABLE 
"RELAY"  OF  THE  RELAY  CONTROL  SLOT 


NO 


YES 


•-' 


CHANGE  THE   BINARY  CODE  TO   "OFF"  VALUE 
SWITCH   =   0 


ASSIGN  THE  NEW  VALUE  TO  THE  VARIABLE 
"RELAY"  OF  THE  RELAY  CONTROL  SLOT 


Figure   3-10.      Flow  chart  of  the  control  program. 
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pumps.   There  were  16  channels  on  the  power  control  board, 
and  the  on/off  status  for  each  channel  was  assigned  by  a  8- 
bit  binary  code  from  the  control  program  executed  in  the 
Zenith  248  microcomputer. 

Final  Control  Elements 

As  stated  in  an  earlier  section,  the  final  control 
elements  were  Cole-Parmer  Master flex  pumps  (Cole-Parmer 
Instrument  Co. ,  Chicago,  Illinois)  to  adjust  the  influent 
and  effluent  flows.   The  pumping  rate  was  set  at  15  ml/min. 
The  required  pumping  time  was  equal  to  the  ratio  of  the 
optimum  flow  volume  divided  by  the  pumping  rate.   In  the 
control  program,  the  resolution  of  the  pumping  time  was  set 
to  "seconds",  and  the  computer  clock  was  used  for  the  real- 
time control . 


CHAPTER  4 
RESULTS  AND  DISCUSSION 


Overall  Performance 


The  operational  characteristics  and  performance  of  the 
experimental  reactor,  including  HRT,  organic  loading  rate, 
methane  content,  gas  production  rate,  methane  production 
rate,  solids,  COD,  volatile  acids,  pH,  organic  removal 
efficiency,  gas  yield  and  methane  yield  are  presented  in 
Figures  4-1,  4-2  and  4-3.   As  shown  in  the  curves  of  the 
methane  content  and  volatile  acids,  instability  was  observed 
at  days  106  and  182.   Nevertheless,  the  reactor  was  able  to 
recover  and  remain  healthy  after  the  feedstock  concentration 
was  changed  and  a  longer  hydraulic  retention  time  was  used. 
This  data  was  helpful  in  setting  the  constraint  for  dynamic 
control,  i.e.,  a  minimum  organic  removal  efficiency  was 
required  to  achieve  the  energy-oriented  objective  functions 
without  sacrificing  the  stability  of  the  reactor.   Based  on 
the  experimental  data,  a  constraint  of  minimum  85%  COD 
removal  efficiency  was  then  selected  and  added  to  the 
objective  functions  for  the  dynamic  control  experiments. 
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Mass  Balance 

To  test  the  reliability  of  the  experimental  data,  mass 
balance  was  conducted  by  using  COD  as  the  index.   The  COD 
mass  flow  of  the  reactor  system  in  a  small  time  period  from 
time  t|(.1  to  time  tk  is  shown  in  Figure  4-4.   From  the 
conservation  law,  mass  balance  for  this  time  block  is 

COD,- (k-1)  +  CODr(k-l)  =  CODr(k)  +  CODe(k-l)  +  CODg(k)   (4-1) 

where  the  subscripts  "i",  "e",  "r"  and  "g"  denote  the  COD 

amount  in  the  influent,  effluent,  reactor  and  gas, 

respectively.   The  corresponding  COD  amount  in  the  gas 

produced  was  calculated  as  follows: 

1  mole     273 

Gas  COD  (g)  =  total  gas  produced  (L)  x  x  

22.4  L      300 

16  g  CH4  44  g  C02 

x  ( x  %  CH4  +  x  %  C02)      (4-2) 

1  mole  CH4  1  mole  C02 

where  the  ratio  273/300  is  the  adjustment  for  the  gas  volume 
since  the  measurement  was  conducted  at  a  room  temperature  of 
27  °C  (=  300  °K) .  Using  eguation  4-1,  the  mass  balance  from 
time  t0  to  tn  was  deduced  as  follows: 

CODj(O)    +  CODr(0)    =  CODr(l)    +  CODJO)    +  COD  (1) 
COD,.  (1)    +  CODr(l)    =  CODr(2)    +  CODe(l)    +  CODg(2) 


COD,- (n-2)  +  CODr(n-2)  =  CODp(n-l)  +  CODe(n-2)  +  COD  (n-1) 
+  COD,- (n-1)  +  CODr(n-l)  =  CODr(n)    +  CODe(n-l)  +  CODg(n) 

EoCODj(t)  +  CODp(0)    =  C0Dr(n)    +  2  C0De(t)  +  2  CODe(t) 

(4-3) 
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70 
or  could  be  stated  as 


Total  input  COD  +  Initial  COD  in  reactor  =  Final  COD  in 
reactor  +  Total  output  COD  +  Total  Gas  COD 

(4-4) 


The  results  listed  in  Table  4-1  show  that  the  ratios  of 
input/output  are  within  an  range  from  1.070  to  1.086  for 
each  particular  testing  period  and  average  from  1.045  to 
1.066  for  the  daily  mass  balance. 

COD-TOC-VS  Correlation 

Generally,  for  a  given  substrate,  a  relationship  exits 
between  the  indexes  of  organic  content  such  as  chemical 
oxygen  demand  (COD) ,  volatile  solids  (VS)  and  organic  carbon 
(TOC) .   In  this  study,  TOC  was  selected  as  the  monitored 
variable  during  dynamic  control  operation  because  it  can  be 
determined  much  faster  than  COD  and  VS.   However,  since 
parameter  values  were  not  available  for  a  TOC-based  dynamic 
model,  it  was  necessary  to  choose  either  COD  or  VS  as  the 
indirect  measurement  which  would  be  determined  by 
correlation  with  the  TOC  measurements.   To  determine  which 
one  was  more  suitable  for  the  secondary  measurement,  a 
comparison  was  made  using  the  experimental  data  in  phase  1. 
The  linear  regression  eguations  were  found  as  follows 
(Figures  4-5  and  4-6) : 

COD  (mg/L)  =  3.15  x  TOC  (mg/L)  +  590,  r2  =  0.99 

(4-5) 

and   VS  (mg/L)   =  2.23  x  TOC  (mg/L)  +  640,  r2  =  0.89 

(4-6) 
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The  result  indicated  that  COD  and  TOC  had  a  high  correlation 
with  an  r  square  of  0.99,  while  the  correlation  between  VS 
and  TOC  was  lower  and  had  an  r  square  of  0.89.   Therefore, 
COD  was  chosen  as  the  state  variable  to  represent  the 
organic  level  for  monitoring  the  reactor  during  dynamic 
control  operation. 

Steady  State  Operation 

During  steady  state  operation  of  the  reactor,  data  were 
collected  primarily  for  use  in  verifying  the  parameter 
estimator  and  constructing  the  initial  model  for  dynamic 
control  operation.   After  recovering  from  a  near- failure 
state,  the  reactor  was  operated  at  different  loading  rates 
with  a  constant  concentration  of  feedstock  and  stayed 
healthy  for  the  remainder  of  the  experiment.   The  loading 
rates  were  applied  in  an  ascending  order  by  gradually 
increasing  the  flow  rate  so  that  shock  loading  of  the 
reactor  would  not  occur.   These  tests  were  described  in  the 
previous  chapter  as  S1#  S2  and  S3.   The  operational 
performance  of  these  tests  are  presented  in  Table  4-2  and  in 
Figures  4-7,  4-8  and  4-9.   These  data  were  helpful  in 
choosing  the  set  point  for  the  objective  function  in 
dynamic  control.   As  shown  in  Table  4-2,  among  the  three 
steady  state  experiments,  S3  was  operated  at  an  average  HRT 
of  21.6  days  with  the  highest  MPR  of  0.487  L  CH4/L-day  and  a 
methane  yield  of  0.317  L  CH4/g  COD  added  (or  0.369  L  CH4/g 
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Table  4-2.   Performance  data  for  steady  state  operation. 


Test 


Parameters 


•1 


Q,  L/day 

HRT,  days 

LRC, 
g  cellulose/L-day 

LRO,    g   COD/L-day 

PH 

COD,    g/L 

TOC,    g/L 

TS,    g/L 

VS,    g/L 

VA,    g/L 

Methane  content,  % 

GPR,  L  gas/L-day 

MPR,  L  CH4/L-day 

YGCb,  L  gas/ 
g  cellulose  added 

YMCb,  L  CJLy 
g  cellulose  added 

YGOb,  L  gas/ 
g  COD  added 

YMOb,  L  CH4/ 
g  COD  added 

COD  removal  eff . ,  % 

TS  removal  eff. ,  % 

VS  removal  eff.,  % 


0.258±0.006a  0.383±0.007   0.463±0.013 
38.8±0.9      26.1±0.5      21.6±0.6 
0.580±0.019   0.862±0.017   1.041±0.029 


0.780±0.019 

7.23±0.07 

3.60±0.30 

0.94±0.10 

7.46±0.27 

2.86±0.25 

0.11±0.03 

51.9±1.3 

0.488±0.053 

0.25410.031 

0.76610.087 


1.16010.022 

7.1710.08 

3.1310.37 

0.8710.12 

7.5510.39 

3.1410.33 

0.1110.04 

51.811.1 

0.76810.039 

0.39810.021 

0.80410.056 


1.40010.039 

7.1410.04 

4.2810.67 

1.1910.23 

8.2510.41 

3.7610.52 

0.2110.12 

51.611.1 

0.94410.086 

0.48710.044 

0.82710.076 


0.39910.051   0.41710.030   0.42610.038 


0.57010.065   0.59810.042   0.61510.056 


0.29710.037   0.30910.023   0.31710.028 


88.111.0 
73.211.0 
87.611.1 


89.611.2 
72.911.4 
86.411.4 


85.912.2 
70.411.5 
83.712.2 


a.  average  1  standard  deviation. 

b.  Values  reported  at  standard  temperature  (0  °C) . 
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79 
COD  destroyed) ,  which  was  close  to  the  theoretical  value  of 
0.35  L  CH4/g  COD  destroyed  reported  by  McCarty  (1964a).   In 
addition,  as  stated  previously,  the  reactor  had  experienced 
instability  at  HRT  of  20  days.   Based  on  this  information, 
the  maximum  MPR  that  could  be  achieved  without  causing 
failure  of  the  reactor  should  be  close  to  the  MPR  of  S3,  and 
a  set  point  higher  than  that  may  cause  instability  of  the 
reactor.   Thus,  a  MPR  of  0.4  L  CH4/L-day,  which  was  about 
80%  of  the  theoretical  maximum  MPR,  was  selected  as  the  set 
point  to  be  tested.   The  COD  removal  efficiency  was  also 
incorporated  to  regulate  the  operation  of  the  reactor. 
Although  the  COD  removal  efficiency  was  the  highest  in  S2  at 
89.6%  and  lowest  in  S3  at  85.9%,  all  three  tests  were  stable 
during  operation.   Therefore,  the  reactor  should  have  been 
capable  of  maintaining  its  stability  if  the  removal 
efficiency  was  85%  or  higher.   This  condition  was  then  added 
to  the  objective  function  as  a  constraint  in  dynamic 
control . 

Verification  of  the  Parameter  Estimator 

Comparison  between  the  Observed  and  the  Predicted  Results 
To  make  sure  that  the  parameter  estimator  functioned 
properly  in  the  dynamic  control  stage,  verification  was 
conducted  by  comparing  the  observed  and  simulated  results. 
A  set  of  parameter  values,  as  listed  in  Table  4-3,  was  found 
using  the  algorithm  and  program  described  in  Chapter  3  and 
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Table  4-3.  Parameter  values  found  by  using  the 

parameter  estimator  and  the  experimental 
data  from  S,  and  S, . 


Parameters 

Values 

Ya,    g  Xa/g 

s 

0.1679 

Yb,    g  Va/g 

s 

0.0703 

Yc,    g  Xm/g 

Va 

0.1815 

Yd,    L  CH4/g 

Va 

2.623 

Ye,    L  C02/g 

Xa 

0.6069 

Yf,    L  C02/g 

Xm 

7.370 

Ks1f    g/L 

2.596 

Ks2,    g/L 

0.7553 

Kdv    day"1 

0.1880 

Kd2,    day"1 

0.03500 

Mm1,    day"1 

0.5361 

Mnc,    day"1 

1.114 

Xa0,    g/L 

1.123 

Xm0,    g/L 

0.8951 

81 
the  experimental  data  from  S2  and  S3.   With  these 
parameters,  the  simulation  results  were  generated  and 
compared  with  the  experimental  data  of  S2  and  S3  in  the 
first  phase  (Figures  4-10  and  4-11) .   The  comparisons  of  the 
average  performance  for  different  testing  periods  are  listed 
in  Table  4-4.   The  results  showed  that  the  curve  of  the 
simulated  data  can  follow  the  trend  of  the  experimental 
data.   This  indicated  that  the  reactor  performance  could  be 
well  predicted  by  using  the  parameter  estimator. 

Different  Numbers  of  Runs  for  the  Parameter  Estimator 

There  was  a  guestion  as  to  whether  the  initial  feasible 
guess  would  affect  the  efficiency  of  the  parameter  estimator 
in  reaching  the  "real"  (or  optimum)  parameter  set.   Also,  if 
the  previously  found  parameter  set  was  used  as  the  initial 
guess  to  start  the  searching,  the  newly  found  parameter  set 
might  fit  the  data  better  than  the  previous  one.   Therefore, 
for  a  group  of  given  data,  more  repetitions  of  the  above 
procedure  would  produce  more  accurate  results,  i.e.,  a 
smaller  residual  sum  of  sguares  (RSS) .   However,  the 
complicated  iteration  procedures  reguired  a  large  amount  of 
computer  time  to  execute  the  program.   Seventy-five  data 
points  and  166  days  of  simulation  time  took  around  3  hours 
of  CPU  time  to  complete  one  run.   Considering  the  sampling 
period  of  12  hours  and  the  time  reguired  to  execute  the 
optimizer,  three  runs  was  the  maximum  that  could  be  adopted. 
A  comparison  of  the  parameters,  total  residual  sum  of 
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85 
squares  and  simulated  results  using  different  numbers  of 
runs  for  the  parameter  estimator  are  displayed  in  Tables  4-5 
and  4-6  and  in  Figures  4-12  and  4-13.   Although  the  tenth 
run  had  the  smallest  RSS  and  was  the  closest  to  the 
experimental  data,  the  difference  was  insignificant  and 
could  be  neglected.   Therefore,  parameters  found  in  the 
third  run  of  the  parameter  estimator  were  used  for  executing 
the  optimizer. 

Dynamic  Control  Using  NSTR 

To  start  a  control  action,  the  initial  model  used  a 
database  which  included  only  the  data  points  in  test  S3  to 
save  computer  time,  and  then  a  new  data  point  was  added  to 
the  database  every  12  hours. 

Since  the  cellulose  is  insoluble  and  degrades  slowly 
compared  with  other  organic  materials  such  as  glucose  (Noike 
et  al.,  1985;  Lee  and  Donaldson,  1984),  there  may  have  been 
a  lag  in  reactor  response  to  the  manipulated  input.   To 
examine  the  proper  response  time  for  this  system,  the 
optimizer  was  initially  designed  to  find  an  optimal  flow 
rate  with  which  the  reactor  MPR  would  reach  the  set  point 
MPR  within  the  next  60  hours.   Different  target  times,  36, 
24  and  12  hours,  were  also  tested,  and  the  output  data  are 
presented  in  Table  4-7  and  in  Figures  4-14  to  4-17. 

The  results  showed  that  the  outputs  for  12  and  60  hours 
were  within  5%  of  the  set  point  after  the  target  time.   In 
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Table  4-5.   Parameter  values  and  RSS  found  by  using 
the  different  numbers  of  runs  for  the 
parameter  estimator. 


First 

Third 

Tenth 

run 

run 

run 

Parameters : 

Ya, 

g  Xa/g 

s 

0.1683 

0.1679 

0.1704 

Yb, 

g  Va/g 

s 

0.0699 

0.0703 

0.0695 

Yc, 

g  Xm/g 

Va 

0.1840 

0.1815 

0.1735 

Yd, 

L  CH4/g 

Va 

2.622 

2.623 

2.676 

Ye, 

L  C02/g 

Xa 

0.6096 

0.6069 

0.6193 

Yf, 

L  C02/g 

Xm 

7.372 

7.370 

7.286 

Ks1f 

g/L 

2.605 

2.596 

2.701 

Ks2, 

g/L 

0.7494 

0.7553 

0.7965 

Kd,, 

day"1 

0.1884 

0.1880 

0.1891 

Kd2, 

day"1 

0.03505 

0.03500 

0.03471 

Hmy 

day"1 

0.5385 

0.5361 

0.5486 

Mm2' 

day"1 

1.113 

1.114 

1.110 

Xa0, 

g/L 

1.127 

1.123 

1.100 

Xm0, 

g/L 

0.8950 

0.8951 

0.8839 

RSS  2.807      2.753      2.680 


87 


P 
C 

o> 
u 

Q) 
4-1 
4-1 
•H 
X) 

C 
<D 
0) 

t 

Q)     • 

O  >H 
C  O 
rd  -P 

grd 
S 

O  -H 
4-1  -p 

u  w 

0)  a) 

a 

T(  Q) 
O  -P 

■p  0) 
(0  6 


(0 

u 
id 


0) 

a)  x! 
cp-p 

(0 

M  4-1 
CD  o 

> 

(0  en 
c 

4-1  3 
O   >h 

C  4-1 

o  o 

W 

•H  W 
*H   ^ 

«s  <u 
ax? 
e  1 
o  B 

u  c 


l 


a) 

H 

id 


Q 
O 

u 


rd 

I- 


0) 


0) 


Q  J 


u 


en 


X!  Z  <* 

S  8 


«  X  <d 


pc;  (d  id 
o       I 


o 


\D 


O 

• 

O 
ID 


H 

rH 


in 


rH 

o 

CM 


■^       in 
o      o 


CM 


n 


en 

en 

en 

r^ 

r- 

r^ 

CO 

CO 

CO 

CO 

CO 

CO 

vo 


n 


>£i 


en 


n 


en 


en 


CO 


in 

• 

o 
in 


CO 


n 


CO 
CO 


n 


en 


en 


CO 


r> 


CM 

o 
in 


H 

n 

IT) 

en 

o 

m 

O 

o 

O 

^r 

LD 

in 

^r 

•** 

^f 

^ 

^ 

"^ 

en 

t-~ 

CM 

o 

r-~ 

CM 

o 

o 

o 

H 

o 

O 

CO 

CO 

CO 

en 

en 

en 

• 

• 

• 

• 

• 

■ 

o 

o 

o 

o 

o 

o 

^r 

CO 

■<* 

H 

n 

^f 

I 

c 

c 

C 

| 

c 

C 

C 

3 

3 

3 

p 

3 

p 

CM 

u 

u 

U 

in 

U 

>H 

u 

CM 

p 

T5 

£1 

CO 

■p 

TS 

si 

to 

U 

-P 

CO 

h 

p 

>1 

^ 

•H 

C 

>1 

CO 

n 

-H 

B 

-H 

-C 

CD 

-H 

£1 

0) 

Cn 

En 

En 

pK 

H 

Eh 

— * 

y— ' 

(\J 

1*1 

CO 

CO 

88 


CO 
(A 


£ 


eg 
6 


o 

1/6 

sppe  e|!je|OA 


o 
to 

I 

c  — 

Q 
-P  O 
C  U 

a>  a 

t4 


IH 

•H 
TS 

O 


Q 
O 
U 


(0  u 

-P  o 

H  -P 

W  g 

Q)  -H 

^  -P 
(0 

73  0) 

a) 

-P  M 

(0  <D 


-P 

0) 

g 

(0 

<t-i  a 


3 

•H 


C  £ 

O  P 

n 

^  o 
(0 

a  in 

e  c 

o  3 

U  V4 


CM 

H 

I 

0) 

3 
tn 

•H 


■JJ9  |BA01U8J  000 


89 


v> 


<M 

O 

03 

u 


P  u 

C  0< 
0) 

u    - 
a)  a 

•H 

T3 


U 

o 

03 
P 


Oh 

U 


M 

o 
-p 

(0 


3 
03 
0 
M  P 
03 
-C   d) 

a) 
■P  fc 

(0  <D 


p 
0) 

e 
u 

(0 


g 

B 

03 


%  Aep-n/T  Aep-n/n 

luo}uoo  eueiue|/\|  ajej  pojd  eueiuew    ©iej  pojd  seg 


0) 

C  £ 

O  P 
03 

•H  4-1 

M  O 
(0 

a  03 

e  c 

o  3 

u  ^ 


m 


(1) 
■H 


90 


Table  4-7.  Controlled  outputs  for  operating  at  different 
target  times  with  a  set  point  of  0.4  L  CH4/L- 
day. 


Output 
variable 

Target  t 

ime 

,    hours 

MPR 
L  CH4/L-day 

60 

36 

24 

12 

Average 
std.  dev. 

0.409 
0.007 

0.423 
0.035 

0.401 
0.016 

0.400 
0.015 

addition,  the  output  for  12  hours  was  closest  to  the  set 
point.   Therefore,  the  reactor  could  be  operated  at  a  target 
time  of  12  hours,  i.e.,  no  lag  of  the  reactor  response  to 
the  manipulated  input.   This  was  coincident  with  the 
original  concept  of  the  control  algorithm,  i.e.,  applying 
the  currently  found  optimum  flow  rate  would  produce  a 
desired  MPR  by  the  next  sampling  time.   The  operational 
criterion  of  12  hours  target  time  was  then  used  throughout 
the  dynamic  control  operation. 

As  described  in  Chapter  3,  there  were  3  tests,  D1#  D2 
and  D3,  in  the  dynamic  control  stage.   The  operational 
objectives  for  these  tests  were:  set  point  control  of  0.4 
CH4/L-day  for  D1 ,  maximum  MPR  operation  for  D2,  and  a  set 
point  of  0.2  CH4/L-day  for  D3.    The  performance  data  for 
these  three  tests  are  presented  in  Table  4-8  and  in  Figures 
4-18,  4-19  and  4-20. 
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Table  4-8.   Performance  data  for  dynamic  control  operation. 


Test 


Parameters 


D 


D, 


Q,  L/day 

HRT ,  days 

LRC, 
g  cellulose/L-day 

LRO,    g  COD/L-day 

PH 

COD,    g/L 

TOC,    g/L 

TS,    g/L 

VS,    g/L 

VA,    g/L 

Methane  content,  % 

GPR,  L  gas/L-day 

MPR,  L  CH4/L-day 

YGCb,  L  gas/ 
g  cellulose  added 

YMCb,  L  CH4/ 
g  cellulose  added 

YGOb,  L  gas/ 
g  COD  added 

YMOb,  L  CH4/ 
g  COD  added 

COD  removal  eff . ,  % 

TS  removal  eff. ,  % 

VS  removal  eff.,  % 


0.398±0.017a  0.452±0.071   0.200±0.019 
25.2±1.1      22.9±4.8      50.5±5.6 
0.895±0.039   1.016±0.161   0.451±0.042 


1.204±0.052 

7.18±0.02 

4.08±0.05 

1.51±0.04 

7.92±0.10 

4.12±0.17 

0.18±0.01 

51.0±0.7 

0.784±0.029 

0.400±0.015 

0.797±0.019 


1.367±0.216 

7.17±0.04 

3.70±0.42 

1.23±0.19 

7.59±0.19 

3.95±0.24 

0.14±0.02 

52.0±1.1 

0.946±0.077 

0.492±0.042 

0.865±0.135 


0.606±0.056 

7.37±0.05 

2.98±0.08 

1.02±0.03 

7.19±0.10 

3.25±0.19 

0.1110.01 

53.0±0.8 

0.394±0.055 

0.209±0.028 

0.804±0.134 


0.407±0.010   0.450±0.070   0.426±0.068 


0.592±0.015   0.643±0.100   0.598±0.099 


0.302±0.007   0.334±0.052   0.317±0.051 


86.5±0.2 
71.6±0.3 
82.210.8 


87.811.4 
72.810.7 
82.911.0 


90.210.3 
74.210.4 
86.010.8 


a.  average  1  standard  deviation. 

b.  Values  reported  at  standard  temperature  (0  °C) . 
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Designated  MPR  Operation 

As  shown  in  Table  4-8,  the  average  MPR  of  test  D,  was 
0.400  L  CH4/L-day  which  was  the  same  as  the  desired  set 
point  and  a  standard  deviation  of  0.015  showed  that  the 
fluctuation  was  very  small.   In  addition,  86.5%  COD  removal 
efficiency,  0.18  g/L  of  volatile  acids  and  a  pH  value  of 
7.18  indicated  that  the  reactor  was  quite  stable  and  well 
controlled. 

The  average  MPR  for  test  D3  was  0.209  L  CH4/L-day  which 
was  nearly  equal  to  the  desired  set  point  0.2  L  CH4/L-day. 
The  reactor  also  remained  stable  with  the  COD  removal 
efficiency,  volatile  acids  and  pH  at  90.2%,  0.11  g/L  and 
7.37,  respectively. 

Maximum  MPR  Operation 

A  major  premise  for  any  operational  objective  is  to 
operate  the  reactor  safely  and  not  to  cause  any  accumulation 
of  toxic  substances.   Therefore,  the  maximum  methane 
production  rate  operation  was  conducted  under  a  constraint 
of  COD  removal,  i.e.,  to  approach  the  highest  MPR  within  a 
predetermined  level  of  COD  removal  efficiency.   As  shown  in 
Table  4-8,  the  average  MPR  of  test  D2  was  0.492  L  CH4/L-day, 
which  was  even  higher  than  the  levels  achieved  earlier  under 
steady  state  operation  (S1f  S2  and  S3)  .   Additionally,  a 
fairly  high  COD  removal  efficiency  of  87.8%  and  a  low 
volatile  acids  of  0.14  g/L  indicated  that  the  reactor  was 
stable,  even  when  operated  at  the  "maximum"  MPR. 
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The  above  results  showed  that  this  system  could  operate 
at  a  high  MPR  and  prevent  instability.   When  the  daily 
performance  was  examined  in  more  detail  (refer  to  Figure  4- 
21) ,  it  was  found  that  the  flow  rate  was  raised  to  a  peak  of 
0.646  L/day  (or  HRT  =15.5  days)  at  day  447.5  (the  second 
sampling  time  for  day  447)  and  resulted  in  the  highest  MPR 
of  0.606  L  CH4/L-day  at  the  next  sampling  time  (day  448). 
In  response  to  this  situation,  the  system  predicted  that  the 
flow  rate  should  be  decreased  to  0.244  L/day  (or  HRT  =  41.0 
days) ,  which  was  the  lowest  point  during  the  maximum  MPR 
operation,  at  the  next  feeding  time  to  prevent  a  decrease  in 
COD  removal  efficiency.   This  helps  explain  why  the  system 
could  remain  stable  under  maximum  MPR  operation. 

Adaptability  of  Changing  Operations 

Besides  its  effectiveness  in  set  point  and  maximum  MPR 
control,  the  control  system  exhibited  excellent 
characteristics  when  changing  operational  objectives.   As 
depicted  in  the  transient  behavior  form  D2  to  D3  (Figure  4- 
22),  it  took  only  3  sampling  periods  (or  1.5  days)  to  bring 
the  MPR  from  a  level  of  around  0.5  to  around  0.2  L  CH4/L- 
day,  and  the  reactor  remained  stable.   After  the  transient 
period,  the  reactor  was  able  to  be  controlled  at  a  set  point 
of  0.2  L  CH4/L-day  of  MPR.   Although  the  rate  of  change  of 
the  output  variable  was  high  at  a  value  of  0.165  L  CH4/L-day 
(or  34%)  per  day,  the  COD,  volatile  acids  and  pH  remained  at 
stable  levels. 
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Gas  Production  Pattern 

The  distribution  of  gas  production  for  two  different 
sampling  periods,  day  427  -  427.5  and  day  440  -  440.5,  was 
observed,  and  the  results  are  shown  in  Figures  4-23  and  4- 
24,  respectively.   As  shown  in  Figure  4-23,  gas  was  produced 
slowly  in  the  first  hour  and  then  increased  steeply  in  the 
second  hour  to  the  highest  rate  of  0.52  L/hr.   Following  the 
peak  production,  the  rate  gradually  decreased  till  the  end 
of  this  period,  and  then  it  started  another  cycle. 
Similarly,  gas  production  for  another  sampling  period 
(Figure  4-24)  increased  slowly  at  the  beginning  and  reached 
the  highest  production  rate  of  0.46  and  0.24  L/L-day  for  gas 
and  methane,  respectively,  in  the  second  hour,  and  then 
gradually  decreased.   Since  gas  production  was  related  to 
the  utilization  of  the  substrate,  these  distributions  curves 
implied  that  a  short  sampling  period  might  help  to  reduce 
the  effect  of  intermittent  feeding.   Therefore,  if  the  high- 
speed computer  and  rapid  response  sensors  are  available, 
feeding  and  sampling  freguency  should  increase  to  improve 
the  smoothness  of  the  operation. 
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CHAPTER  5 
SUMMARY  AND  CONCLUSIONS 


A  computerized  system  using  a  nonlinear  self-tuning 
regulator  (NSTR)  was  developed  for  the  operation  of  an 
anaerobic  continuously  stirred  tank  reactor.   The  NSTR 
control  algorithm  was  the  combination  of  an  adaptive 
parameter  estimator  and  an  optimizer.   Based  on  historical 
data,  the  parameter  estimator  found  a  set  of  "best-fitted" 
parameters,  which  were  then  used  by  the  optimizer  to  locate 
the  optimum  manipulated  input  (flow  rate)  to  accomplish  the 
desired  output  for  the  next  feeding. 

A  two-culture  dynamic  model,  implemented  in  the  NSTR 
control  algorithm,  was  developed  to  describe  the  microbial 
interactions  of  an  anaerobic  continuously  stirred  tank 
reactor.   This  model  was  able  to  predict  the  dynamic 
behavior  of  the  reactor,  including  concentrations  of  organic 
substrate,  volatile  acids,  acetogens  and  methanogens,  and 
also  the  gas  production  rate  and  gas  composition. 

To  verify  the  proposed  control  algorithm,  a  bench-scale 
reactor  eguipped  with  the  computer  control  hardware  was  set 
up  and  operated  for  different  objectives.   Results  of  the 
controlled  output  showed  that  the  control  algorithm  could  be 
successfully  employed  in  this  computer  control  system.   In 
one  of  the  set  points,  at  a  methane  production  rate  (MPR)  of 
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0.4  L  CH4/L-day,  the  average  output  was  the  same  as  the 
desired  value  with  very  little  fluctuation.   For  another  set 
point  control,  the  average  output  of  0.209  L  CH4/L-day  was 
only  slightly  higher  than  the  desired  value  of  0.2  CH4/L- 
day.   Also,  in  control  for  maximum  MPR  operation,  the 
reactor  was  able  to  perform  at  a  high  methane  production 
rate  and  remain  stable.   The  COD  removal  efficiency  for  all 
three  tests  was  higher  than  the  predetermined  minimum  limit 
of  85%. 

Besides  its  capability  in  controlling  the  methane 
production  rate,  this  system  also  showed  an  excellent 
adaptability  in  changing  the  operational  objectives.   It 
took  only  3  sampling  periods  (1.5  days)  to  change  the  MPR 
from  a  level  of  around  0.5  to  0.2  L  CH4/L-day,  and  the 
reactor  remained  stable. 

In  general,  the  advantages  of  a  NSTR  control  algorithm 
can  be  summarized  as  follows: 

1.  Although  an  extremely  detailed  model  is  good  in  providing 
more  information  (without  considering  computer  execution 
time) ,  it  is  not  a  reguisite  for  this  control  system, 
since  the  parameter  estimator  will  search  the  optimum 
parameter  set  based  on  the  observed  data. 

2.  If  there  are  not  enough  data  to  generate  a  good-fitting 
parameter  set,  it  will  gradually  approach  the  real 
process  by  adapting  from  the  increasing  amount  of  data. 

3.  Although  the  NSTR  in  this  study  was  designed  especially 
for  anaerobic  processes,  a  NSTR  control  algorithm  is  also 
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applicable  to  wherever  the  dynamic  model  for  that 
particular  process  is  available. 
The  following  conclusions  can  be  drawn  from  this  study: 

1.  A  two-culture  dynamic  model  was  developed  which  was  able 
to  predict  the  dynamic  behavior  of  the  system. 

2 .  A  NSTR  control  algorithm  was  developed  which  included  a 
parameter  estimator  and  an  optimizer.   The  optimizer  was 
able  to  predict  the  optimum  flow  rate  to  achieve  the 
desired  output. 

3.  A  bench-scale  reactor  was  successfully  set  up  and 
operated  for  verifying  the  proposed  control  algorithm. 

4.  At  a  set  point  of  0.4  L  CH4/L-day,  the  average  output  was 
the  same  as  the  desired  value  with  very  small 
fluctuation. 

5.  At  a  set  point  of  0.2  L  CH4/L-day,  the  average  output  was 
only  slightly  higher  than  the  desired  value. 

6.  At  the  maximum  MPR  operation,  the  average  output  was  the 
highest  in  all  of  the  tests  and  remained  stable. 

7.  The  COD  removal  efficiency  for  all  three  tests  in  dynamic 
control  was  higher  than  the  predetermined  minimum  limit 
of  85%. 


CHAPTER  6 
RECOMMENDATIONS  FOR  FURTHER  STUDY 


The  results  from  this  study  are  promising.   However, 
there  is  still  room  for  improvement  in  the  control  of 
anaerobic  reactors.   The  following  recommendations  are  made 
for  further  research: 

1.  To  test  the  NSTR  control  system  in  a  pilot  scale  reactor. 

2 .  To  test  the  NSTR  control  system  at  variable  feeding 
conditions  such  as  different  types  of  substrate  and 
variable  feedstock  concentrations. 

3 .  To  develop  a  control  system  based  on  the  NSTR  which  is 
applicable  to  other  types  of  substrates  and  reactors. 

4 .  When  sensors  and  fast-determination  methods  are 
available,  include  more  variables  for  calculating  the 
total  residual  sum  of  squares  in  the  parameter  estimator. 

5.  To  develop  a  more  accurate  process  model  including  more 
microbial  interactions  such  as  hydrolysis  and 
homoacetogenic  steps.   This  will  shorten  the  convergence 
time  when  estimating  parameters  and  require  less  data  for 
reaching  the  real  process. 

6.  To  develop  a  model  in  which  the  index  of  organic  level 
can  be  measured  on-line,  such  as  TOC. 

7.  To  increase  the  feeding  and  sampling  frequency  for 
improving  the  smoothness  of  control. 
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8.  To  cut  the  computer  time  by  improving  the  efficiency  of 
the  parameter  estimator  and  optimizer. 

9.  To  investigate  the  possibility  of  incorporating  expert 
system  technigues  in  the  control  system  to  increase  its 
flexibility,  such  as  dealing  with  multiple  objectives  or 
multiple  constraints. 


APPENDIX  A 
NOMENCLATURE 


COD  Chemical  oxygen  demand  concentration  in  reactor, 
g/L 

COD0  Initial  COD  concentration,  g/L 

COD,-  Influent  COD  concentration,  g/L 

GPR  Volumetric  biogas  production  rate,  L/L-day 

HRT  Hydraulic  retention  time,  days 

Kd1  Decay  rate  coeff.  for  Xa,  day* 

Kd2  Decay  rate  coeff.  for  Xm,  day" 

Ks1  Saturation  coeff.  for  Xa,  g/L 

Ks2  Saturation  coeff.  for  Xm,  g/L 

LRC  Cellulose  loading  rate,  g  cellulose/L-day 

LRO  COD  loading  rate,  g  COD/L-day 

MPR  Volumetric  methane  production  rate,  L  CH4/L-day 

MPR  Predicted  methane  production  rate,  L  CH^/L-day 

MPR  Set  point  of  volumetric  methane  production  rate, 
L  CH4/L-day 

PCH4  Methane  content,  % 

Q  Flow  rate,  L/day 

QCH4  Daily  methane  production,  L  CH4/day 

Qco2  Daily  carbon  dioxide  production,  L  C02/day 

QGAS  Daily  biogas  production,  L/day 

Qm  Manipulated  flow  rate,  L/day 

S  Substrate  concentration  in  reactor,  g/L 

111 


112 


S0        Initial  substrate  concentration,  g/L 


Si  Influent  substrate  concentration,  g/L 

V  Liquid  volume  of  reactor,  L 

Va  Volatile  acid  in  reactor,  g/L 

Va0  Initial  volatile  acid  in  reactor,  g/L 

Vai  Influent  volatile  acid,  g/L 

VIN  Volume  of  influent  flow,  L 

Xa  Acetogen  concentration  in  reactor,  g/L 

Xa0  Initial  acetogen  concentration  in  reactor,  g/L 

Xai  Influent  acetogen  concentration,  g/L 

Xm  Methanogen  concentration  in  reactor,  g/L 

Xm0  Initial  methanogen  concentration  in  reactor,  g/L 

Xmi  Influent  methanogen  concentration,  g/L 

Ya  Xa  yield  coeff.,  g  Xa  prod./g  s  used 

Yb  Va  yield  coeff.,  g  Va  prod./g  S  used 

Yc  Xm  yield  coeff.,  g  Xm  prod/g  Va  used 

Yd  CH4  yield  coeff.  ,  L  CH4  prod./g  Va  used 

Ye  C02  yield  coeff.  via  acetogen,  L  C02/g  Xa  prod. 

Yf  C02  yield  coeff.  via  methanogen,  L  C02/g  Xm  prod. 

YGC  Cellulose  gas  yield,  L  gas/g  cellulose  added 

YGO  COD  gas  yield,  L  gas/g  COD  added 

YMC  Cellulose  methane  yield,  L  CH4/g  cellulose  added 

YMO  COD  methane  yield,  L  CH4/g  COD  added 

Hy  Specific  growth  rate  of  Xa,  1/day 

M2  Specific  growth  rate  of  Xm,  1/day 

Mmi  Maximum  specific  growth  rate  of  Xa,  day"1 

Mm2  Maximum  specific  growth  rate  of  Xm,  day"1 


APPENDIX  B 
PROGRAM  LISTING 


Parameter  Estimator 


C  ********************************************************** 

c  *  * 

C  *  Filename:  PE.FOR  * 

C  *  * 

C  *  This  program  is  the  first  part  of  the  Nonlinear  Self-  * 

C  *  Tuning  Regulator.   It  estimates  the  parameters  of  the  * 

C  *  anaerobic  digestion  model,  which  will  be  loaded  to  the  * 

C  *  optimizer  (OP. FOR)  for  finding  the  optimum  manipulated  * 

C  *  variable  -  flow  rate.   The  residual  sum  of  squares  is  * 

C  *  chosen  as  the  criteria  for  finding  the  best  set  of  * 

C  *  parameters.   The  searching  procedures  are  based  on  a  * 

C  *  modified  Box  Complex  Algorithm.   When  stopping  rule  is  * 

C  *  met,  a  set  of  parameter  will  be  determined.   Beginning  * 

C  *  time  and  ending  time  need  be  defined  in  the  DATA  * 

C  *  statement:  TINI  and  TFIN.  * 

C  *  * 

C  *  Subroutine  (or  Function)  No. :  * 

C  *    1.  PARA  * 

C  *    2.  RAN3  * 

C  *    3.  COMPLEX  * 

C  *    4.  CALC  * 

C  *    5.  INIT  * 

C  *    6.  RUNGE  * 

C  *    7.  OUTPUT  * 

C  *    8.  SUMSQR  * 

C  *    9.  RATE  * 

C  *  * 

C  *  Subroutines  called  in  main  program:  PARA,  COMPLEX  * 

C  *  * 

C  *  Substrate  used:  COD  * 

C  *  * 

C  *  Variables  (not  including  those  in  nomenclature) :  * 

C  *    A       =  Current  set  of  parameter  * 

C  *   CENTRO  =  Average  of  parameter  values  exclude  the  * 

C  *             worst  guess  * 

C  *   CH4     =  Predicted  methane  produced  in  one  time  * 

C  *             step,  L  * 

C  *   CH4T    =  Predicted  methane  produced  in  one  sampling  * 

C  *             period,  L  * 

C  *   C02     =  Predicted  carbon  dioxide  produced  in  one  * 

C  *             time  step,  L  * 
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c 

* 

C02T 

c 

* 

c 

* 

CONS I Z 

c 

* 

DELT 

c 

* 

G 

c 

* 

GAS 

c 

* 

c 

* 

GAST 

c 

* 

c 

* 

GP 

c 

* 

c 

* 

IBEGIN 

c 

* 

IC 

c 

* 

c 

* 

IEND 

c 

* 

IR 

c 

* 

c 

* 

ISAMP 

c 

* 

MAXIM 

c 

* 

MI 

c 

* 

MINIM 

c 

* 

MP 

c 

* 

c 

* 

MX1 

c 

* 

MX  2 

c 

* 

NEQN 

c 

* 

NPRAM 

c 

* 

NT 

c 

* 

OUTFILE 

c 

* 

PRDEL 

c 

* 

PQ 

c 

* 

R 

c 

* 

c 

* 

RSS 

c 

* 

c 

* 

SCOD 

c 

* 

SGPR 

c 

* 

SMPR 

c 

* 

SPCH4 

c 

* 

SSCOD 

c 

* 

SSGPR 

c 

* 

SSMPR 

c 

* 

SSPCH4 

c 

* 

SSVA 

c 

* 

SUMNEW 

c 

* 

SUMDIF 

c 

* 

c 

* 

SVA 

c 

* 

SXA 

c 

* 

SXM 

c 

* 

TFIN 

c 

* 

TIME 

c 

* 

TINI 

c 

* 

Y 

Predicted  carbon  dioxide  produced  in  one 
sampling  period,  L 

Size  of  convergence,  square  root  of  SUMDIF 
Step  size  for  integration,  day 
Guesses  of  parameters 

Predicted  biogas  produced  in  one  time 
step  (CH4  +  C02) ,  L 

Predicted  biogas  produced  in  one  sampling 
period,  L 

Observed  gas  produced  in  one  sampling 
period,  L 

Sample  number  to  start  simulation 
The  set  number  of  parameters  with  the  best  * 

guess  (minimum  RSS)  * 

Sample  number  to  end  simulation  * 

The  set  number  of  parameters  with  the  * 

worst  guess  (maximum  RSS)  * 

Sample  number  * 

Upper  bound  of  parameter  * 

Minimum  RSS  * 

Lower  bound  of  parameter  * 

Observed  methane  produced  in  one  * 

sampling  period,  L  * 

Maximum  RSS  * 

Maximum  RSS  after  adding  the  new  guess  * 

Number  of  rate  equations  * 

Number  of  parameters  to  be  estimated  * 

Counter  to  cumulate  time  step  * 

Filename  to  store  output  * 

Printing  period,  day  * 

Flow  rate  of  previous  sample  * 

Values  of  RSS  for  individual  set  of  * 

parameters  * 

Total  residual  sum  of  squares,  sum  of  * 

SSGPR,  SSMPR,  SSCOD  and  SSPCH4  * 

Predicted  COD  concentration  ,  g/L  * 

Predicted  gas  production  rate,  L/L-day  * 
Predicted  methane  production  rate,  L/L-day  * 

Predicted  methane  content,  %  * 

Residual  sum  of  squares  for  COD  * 

Residual  sum  of  squares  for  GPR  * 

Residual  sum  of  squares  for  MPR  * 

Residual  sum  of  squares  for  PCH4  * 

Residual  sum  of  squares  for  VA  * 

Sum  of  parameter  values  for  all  guesses  * 

Sum  of  squares  of  RSS  difference  for  all  * 

guesses  from  the  CENTRO  * 

Predicted  VA  concentration,  g/L  * 

Predicted  acetogen  concentration,  g/L  * 

Predicted  methanogen  concentration,  g/L  * 

Time  to  end  simulation,  day  * 

Simulation  time,  day  * 

Time  to  start  simulation,  day  * 

State  variables,  g/L  * 
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C  *  YO  =  Initial  conditions  of  state  variables,  g/L  * 
C  *  Note:  Variables  Ydot,  PY1,  PY2,  PY3 ,  PYdotl,  PYdot2  * 
C  *  and  PYdot3  were  used  and  defined  in  subroutine  * 
C  *       RUNGE.  * 

C  *  * 

C  ********************************************************** 

C  * 

C  ***  Main  Program 

C  * 

REAL*8  R(16) ,G(16,14) ,A(14) ,MAXIM(14) ,MINIM(14) ,MX1, 
&   MX2 , MI , SUMNEW ( 14 ) , RSS , CENTRO ( 14 ) , SUMDIF , CONSIZ , 
&   MU1 , MU2 , MUHAT1 , MUHAT2 , Ksl , Ks2 , Kdl , Kd2 , SUMDIF 

REAL  MP,MPR 

CHARACTER  0UTFILE*12 

COMMON  /INPUT/NEQN,OUTFILE,NPARM 

COMMON  /Y0/COD0,VA0,XA0,XM0 

COMMON  /SVAR/Y(10) , Ydot (10) ,COD,Va 

COMMON  /RUNKU1/PY1(10) ,PY2(10) ,PY3(10) 

COMMON  /RUNKU2/PYdotl(10) ,PYdot2(10) ,PYdot3(10) 

COMMON  /TIMER/ISAMP , TFIN , TINI , DELT , PRDEL , TIME , NT 

COMMON  /YVALUE/Ya,Yb,Yc,Yd,Ye,Yf 

COMMON  /KVALUE/Ksl,Ks2,Kdl,Kd2 

COMMON  /GROWTH/MUHAT1 , MUHAT2 , MU1 , MU2 

COMMON  /TANK/Q,V,HRT,VIN 

COMMON  / INFLOW/ CODi,Vai,Xai,Xmi 

COMMON  /GAS/CH4T,GAST,MP,GP,GPR,MPR,PCH4,SGPR, 
&  SMPR,SPCH4 

COMMON  /SS/SSGPR , SSMPR , SSCOD , SSVA , SSPCH4 , RSS , 
&  SCOD,SVA,SXA,SXM 

COMMON  /PARM/A,IC 

COMMON  /GUESS/G 

COMMON  /RANGE/MINIM, MAXIM 

COMMON  /COUNTER/IBEGIN,IEND 

DATA  V/10.0/ 

DATA  CODi/30 . 2595/ , Vai/0 . 37875/ , Xai/0 . 0/ , Xmi/0 . 0/ 

DATA  COD0/2.751/,VA0/0.120/ 

DATA  TINI/347./,TFIN/419.5/,IBEGIN/126/,IEND/167/ 

DATA  DELT/ 0.01/ 

DATA  MINIM/. 1, . 05 , . 03 , . 3 , . 5 , 2 . , . 1 , . 15, . 01, . 01 , 
&  .04, .004, .1, .01/ 

DATA  MAXIM/1,2,1,5,2,10,10,10, . 5, . 05 , 5 , 1 . 5, 2 . , 1./ 

NEQN    =  4 

NPARM   =  14 

N       =  NPARM 

IC      =1 

I COUNT  =  1 

OPEN (UNIT=51 , FILE= ' PARM. OPT » ) 

OPEN (UNIT=52 , FILE= • BESTPARM ' ) 

C  * 

C  ***  Repeat  finding  local  optimum  icount  times,  icount  can 

C  ***  be  changed  as  user's  wish.   Here  icount  is  set  as  1. 
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C  * 

DO  1  II=l,ICOUNT 

PRINT  *,  'II  =  ■ ,11 

PRINT  *,  •    ' 

CALL  PARA 

CALL  COMPLEX 

WRITE (52, 2001)  (G (IC, J) , J=1,N) 
2001  FORMAT ( 14 F9. 6) 
1  CONTINUE 

WRITE(51,2001)  (G(IC,J) ,J=1,N) 
CLOSE (51) 
CLOSE (52) 

STOP 
END 

C  ********************************************************** 
c  *  * 

C  *  #1  Para  * 

C  *  * 

C  *  Subroutine  to  initialize  guesses  of  Ya,  Yb,  Yc,  Yd,  Ye,* 
C  *  Yf,  K2,  Ksl,  Ks2,  Kdl,  Kd2 ,  MUHAT1,  MUHAT2 ,  XaO  and  * 
C  *  XmO.  These  guesses  were  grouped  as  a  16  by  14  array.  * 
C  *  * 

C  ********************************************************** 

SUBROUTINE  PARA 

REAL*8  G ( 16 , 14 ) , MAXIM ( 14 ) , MINIM ( 14 ) , A ( 14 ) 

CHARACTER  OUTFILE*12 

COMMON  /INPUT/NEQN,OUTFILE,NPARM 

COMMON  /Y0/COD0,VA0,XA0,XM0 

COMMON  /GUESS/G 

COMMON  /PARM/A,IC 

COMMON  /RANGE/MINIM, MAXIM 

N      =  NPARM 
ISEED  =  97 

OPEN (UNIT=4  5 , FILE= ' PARM . PE • ) 
READ(45,*)  (G(IC,J) ,J=1,N) 
REWIND (45) 
DO  101  J=1,N 
G(1,J)=G(IC,J) 

101  CONTINUE 

DO  102  J=1,N 

DO  103  1=2, N+2 

G(I,J)=(MINIM(J)+RAN3 (ISEED) * (MAXIM (J) -MINIM (J) ) 
&  +G(IC,J))/2 

103     CONTINUE 

102  CONTINUE 

WRITE(45,2101)  ( (G(I,J) ,J=1,N) ,1=1, N+2) 
2101   FORMAT ( 14 F9. 6) 
CLOSE (45) 

RETURN 
END 
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C  ********************************************************** 
C  *  #2  Ran3  * 

C  *  * 

C  *  A  function  of  random  number  generator  for  uniform  * 
C  *  distribution,  using  the  method  developed  by  Donald  E.  * 
C  *  Knuth  (1981).  The  code  is  from  William  H.  Press  * 
C  *  et  al.  ("Numerical  Recipes",  1986)  * 

C  *  * 

C  ********************************************************** 

FUNCTION  RAN3(IDUM) 

PARAMETER  (MBIG=1000000000 ,MSEED=161803398 , 
&  MZ=0,FAC=l.E-9) 

DIMENSION  MA(55) 
DATA  IFF  /0/ 

IF  (IDUM.LT.0.OR.IFF.EQ.0)  THEN 
IFF=1 

MJ=MSEED-IABS (IDUM) 
MJ=MOD(MJ,MBIG) 
MA(55)=MJ 
MK=1 
DO  201  1=1,54 

II=MOD(21*I,55) 

MA(II)=MK 

MK=MJ-MK 

IF  (MK.LT.MZ)  MK  =  MK  +  MBIG 

MJ  =  MA(II) 

201  CONTINUE 

DO  202  K=l,4 
DO  203  1=1,55 

MA(I)  =  MA(I)  -  MA(1+MOD(I+30,55) ) 
IF  (MA(I) .LT.MZ)  MA(I)=MA(I)+MBIG 
203       CONTINUE 

202  CONTINUE 
INEXT   =  0 
INEXTP  =31 
IDUM    =  1 

ENDIF 

INEXT  =  INEXT  +  1 

IF  (INEXT. EQ. 56)  INEXT=1 

INEXTP=INEXTP+1 

IF  (INEXTP. EQ. 56)  INEXTP=1 

MJ  =  MA (INEXT)  -  MA (INEXTP) 

IF  (MJ. LT.MZ)  MJ=MJ+MBIG 

MA (INEXT)  =  MJ 

RAN3       =  MJ  *  FAC 

RETURN 
END 

C  ********************************************************** 
c  *  * 

C  *  #3  Complex  * 

C  *  * 


This  subroutine  uses  the  Box's  Complex  Algorithm  to 
optimize  the  parameters.   (Box,  M.  J.,  1965,  A  new 
method  of  constrained  optimization  and  a  comparison 
with  other  methods.   Computer  Journal,  Vol.  8,  No.l, 
pp. 42-52. ) 

Need  call  subroutine:  CALC 
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* 
* 
* 
* 
* 
* 
* 
* 


c 

C 

c 
c 
c 
c 
c 
c 

c  ********************************************************** 

SUBROUTINE  COMPLEX 

REAL*8  R(16) ,G(16,14) ,A(14) ,H(16,14) ,MAXIM(14) , 
&  MINIM(14) ,SUMNEW(14) ,CENTR0(14) ,RSS,MX1,MX2 ,MI ,MU1 , 
&  MU2 , MUHAT1 , MUHAT2 , Ksl , Ks2 , Kdl , Kd2 , SUMDIF , CONSIZ 
REAL  MP,MPR 
CHARACTER  0UTFILE*12 
COMMON  /INPUT/NEQN,OUTFILE,NPARM 

/Y0/COD0 , VAO , XAO , XMO 

/SVAR/Y(10) ,Ydot(10) ,COD,Va 

/RUNKU1/PY1(10) ,PY2(10) ,PY3(10) 

/RUNKU2/PYdotl(10) ,PYdot2(10) ,PYdot3(10) 

/TIMER/ ISAMP , TFIN , TINI , DELT , PRDEL, TIME , NT 

/ YVALUE/ Ya , Yb , Yc , Yd , Ye , Yf 

/KVALUE/Ksl , Ks2 , Kdl , Kd2 

/GR0WTH/MUHAT1 , MUHAT2 , MU1 , MU2 

/TANK/Q , V , HRT , VIN 

/INFLOW/ CODi , Vai , Xai , Xmi 

/GAS/CH4T , GAST , MP , GP , GPR , MPR , PCH4 , 
SGPR,SMPR,SPCH4 
COMMON  /SS/SSGPR, SSMPR, SSCOD, SSVA, SSPCH4 , RSS , 
i  SCOD,SVA,SXA,SXM 

/PARM/A,IC 

/GUESS/G 

/RANGE/MINIM , MAXIM 

/COUNTER/ I BEGIN, I END 


COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 


COMMON 
COMMON 
COMMON 
COMMON 


DATA  E1,E2/0.0001, . 1/ 


C 

c 
c 
c 
c 


* 

*** 
*** 
*** 

* 


El  is  a  predetermined  number  for  convergence  test, 
E2  is  another  termination  criterion  for  test  the 
absolute  RSS  value  of  the  latest  calculation. 


RSS  =  0. 
N    =  NPARM 

PRINT  *,  'Initial  guesses:' 
PRINT  *,  '    ' 
DO  301  1=1, N+l 
DO  302  J=1,N 
A(J)  =  G(I,J) 
3  02    CONTINUE 
CALL  CALC 
R(I)  =  RSS 
301  CONTINUE 
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C  * 

C  ***  Calculating  an  RSS  value  for  each  set  of  the 

C  ***  estimates,  making  a  total  of  N+l  RSS.   The  criterion 

C  ***  to  choose  an  optimal  set  of  parameter  values  is  based 

C  ***  on  total  residual  sum  of  squares 

C  * 

C  ***  Initialize  the  maximal  and  minimal  RSS  as  R(l) 

C  * 

320  MX1  =  R(l) 
MI   =  R(l) 

C  * 

C  ***  Selecting  the  worst  and  best  RSS  and  designating  their 

C  ***  indexes  as  IR  &  IC. 

C  * 

DO  303  I=1,N+1 

IF(R(I) .GE.MX1)  THEN 
C  ***   Worst  guess 
MX1  =  R(I) 
IR  =  I 
ENDIF 

IF(R(I) .LE.MI)  THEN 
C  ***   Best  guess 
MI  =  R(I) 
IC  =  I 
ENDIF 
303  CONTINUE 

C  * 

C  ***  Minimum  RSS  test 
C  * 

IF(R(IC) .LT.E2)  THEN 

PRINT  *,  'Pass  minimum  RSS  test1 
GOTO  304 
ENDIF 

C  * 

C  ***  Convergence  test 

C  * 

DO  310  1=1, N 
SUMNEW(I)  =  0. 
DO  311  J=1,N+1 

SUMNEW ( I ) =SUMNEW ( I ) +G ( J , I ) 
311      CONTINUE 

CENTRO(I)  =  (SUMNEW(I)-G(IR,I) )  /  N 
A(I)       =  CENTRO(I) 
310    CONTINUE 
CALL  CALC 
RCEN  =  RSS 
SUMDIF  =  0. 
DO  351  1=1, N+l 

SUMDIF  =  SUMDIF  +  (RCEN-R(I) ) **2 
351    CONTINUE 

CONSIZ  =  DSQRT (SUMDIF/ (N+l) ) 
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IF  (CONSIZ  .GE.  El)  THEN 
C  * 

C  ***  The  following  determines  a  new  set  of  guesses 
C  * 

C  ***  G(N+2,I)  is  a  temporary  cell  to  store  the  new  guess 
C  ***  G(N+2,I)=(l+a)*SUMNEW(I)/N  -  ( (1+a) /N+l) *G (IR, I) 
C  ***  Box  recommended  a=1.3,  a=l  is  used  here 
C  * 

DO  313  I=1,N 

G(N+2,I)=2.*SUMNEW(I)/N  -  (2 ./N+l) *G(IR, I) 

C  * 

C  ***  Explicit  constraints  test  for  the  parameters 

C  * 

IF(G(N+2,I) .LE.MINIM(I) )  G (N+2 , I) =MINIM(I) 
IF(G(N+2,I) .GT.MAXIM(I) )  G (N+2 , I) =MAXIM(I) 
A(I)  =  G(N+2,I) 
313    CONTINUE 

CALL  CALC 
R(N+2)  =  RSS 

IF  (R(N+2) .LT.R(IR) )  THEN 
DO  315  NA=1,N 

G(IR,NA)  =  G(N+2,NA) 

315  CONTINUE 
R(IR)  =  R(N+2) 

MX2    =  DMAX1(R(1) ,R(2) ,R(3) ,R(4) ,R(5) ,R(6) , 
&       R(7),R(8),R(9),R(10),R(11),R(12),R(13),R(14)) 
IF  (DABS (MX2-R(IR) ) .GT. 0.0)  THEN 
OPEN (UNIT=4  5 , FILE= ' PARM . PE ' ) 
WRITE (45, 2304)  ((G(I,J) ,J=1,N) ,1=1, N+2) 
2304        FORMAT (14F9. 6) 
CLOSE (45) 
ELSE 
C  * 

C  ***  Repeat  worst  point,  reducing  the  value  by  half 
C  * 

PRINT  *,  'Repeat  worst  point' 
DO  316  K1=1,N+1 
DO  317  L1=1,N 

G(K1,L1)  =  (G(K1,L1)+G(IC,L1) )  /  2. 
IF(G(K1,L1) .LE.MINIM(Ll) )  G(K1, LI) =MINIM(L1) 
IF(G(K1,L1) .GT.MAXIM(Ll) )  G(K1, LI) =MAXIM(L1) 
A(L1)     -  G(K1,L1) 
317  CONTINUE 

CALL  CALC 
R(K1)  =  RSS 

316  CONTINUE 

OPEN (UNIT=45 , FILE=  »  PARM . PE ' ) 
WRITE (45, 23  04)  ((G(I,J) ,J=1,N) ,1=1, N+2) 
CLOSE(45) 
ENDIF 
ELSE 
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C  * 

C  ***  Reducing  the  values  by  half  before  a  new  flipping 

C  * 

PRINT  *,  'Flipping' 
DO  318  K=1,N+1 
DO  319  L=1,N 

G(K,L)  =  (G(K,L)+G(IC,L))  /  2. 
IF(G(K,L) .LE.MINIM(L) )  G (K, L) =MINIM(L) 
IF(G(K,L) .GT.MAXIM(L) )  G (K, L) =MAXIM(L) 
A(L)    =  G(K,L) 
319        CONTINUE 

CALL  CALC 
R(K)  =  RSS 
318      CONTINUE 

OPEN (UNIT=45 , FILE= ' PARM. PE ' ) 
WRITE (45, 2304)  ((G(I,J) ,J=1,N) , I=l,N+2) 
CLOSE(45) 
ENDIF 

GOTO  320 
ELSE 

WRITE(6,2305) 
2305    FORMAT(10X,'   **  cannot  improve  any  further  ***•) 
ENDIF 

C  * 

C  ***  Print  the  results 

C  * 

304  PRINT  *,'IC(best  guess)=',  IC 
DO  321  1=1, N 
A(I)  =  G(IC,I) 
WRITE(6,*)  I,G(IC,I) 
321  CONTINUE 
CALL  CALC 
WRITE (6 , *)  SSGPR, SSMPR, SSCOD, SSPCH4 ,RSS 

RETURN 
END 

C  ********************************************************** 
c  *  * 

C  *  #4  Calc  * 

C  *  * 

C  *  Subroutine  to  conduct  the  dynamic  methane  fermentation  * 
C  *  and  use  4th  order  Runge-Kutta  method  for  integration.  * 
C  *  * 

C  *  Need  call  subroutines:  INIT,  RUNGE,  OUTPUT,  SUMSQR  * 
C  *  * 

C  ********************************************************** 

SUBROUTINE  CALC 

REAL*8   R(16) ,G(16,14) ,A(14) , MAXIM (14) , MINIM (14) , 
&   MX1 , MX2 , MI , SUMNEW ( 14 ) , RSS , 
&   MUl,MU2,MUHATl,MUHAT2,Ksl,Ks2, 
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&   Kdl , Kd2 

REAL  MP,MPR 

CHARACTER  0UTFILE*12 

COMMON  /INPUT/NEQN,OUTFILE,NPARM 

COMMON  /Y0/COD0,VA0,XA0,XM0 

COMMON  /SVAR/Y(10) ,Ydot(10) ,COD,Va 

COMMON  /RUNKU1/PY1(10) ,PY2(10) ,PY3(10) 

COMMON  /RUNKU2/PYdotl(10) ,PYdot2(10) ,PYdot3(10) 

COMMON  /TIMER/ ISAMP , TFIN , TINI , DELT , PRDEL , TIME , NT 

COMMON  /YVALUE/Ya,Yb,Yc,Yd,Ye,Yf 

COMMON  /KVALUE/Ksl,Ks2,Kdl,Kd2 

COMMON  /GROWTH/MUHAT1 , MUHAT2 , MU1 , MU2 

COMMON  /TANK/Q , V , HRT , VIN 

COMMON  /INFLOW/ CODi , Vai, Xai, Xmi 

COMMON  /GAS/CH4T,GAST,MP,GP,GPR,MPR,PCH4, 
&  SGPR,SMPR,SPCH4 

COMMON  /SS/SSGPR, SSMPR, SSCOD , SSVA , SSPCH4 , RSS , 
&  SCOD,SVA,SXA,SXM 

COMMON  /PARM/A,IC 

COMMON  /COUNTER/ I BEGIN, I END 

CALL  INIT 

C  * 

C  ***  Update  state  variables 

C  * 

Y(l)  =  CODO 

Y(2)  =  VAO 

Y(3)  =  XAO 

Y(4)  =  XMO 

OPEN  (UNIT=46,FILE=,CSIMUPE.OUTl) 

C  * 

C  ***  CSIMUPE.OUT  is  the  filename  for  simulated  performance 
C  ***  data.   COPT. DAT  is  the  filename  of  observed  data. 
C  * 

PDAY   =  TINI 
PQ     =  0.452 

OPEN  (UNIT=47,FILE=' COPT. DAT') 
407  Q  =  PQ 

READ(47,1401)  ISAMP, RDAY, PQ, VIN, GP, MP, GPR, 
&   MPR,PCH4,COD,VA 
1401  FORMAT(1X,I3,F6.1,2F6.3,2F6.2,2F6.3,F6.2,2F7.0) 
PRDEL  =  RDAY 

C  * 

C  ***  Adjust  unit  to  g/L 

C  * 

COD    =  COD  /  1000. 

VA     =  VA  /  1000. 

C  * 

C  ***  Finding  the  concentration  of  state  variables  at  the 
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C  ***  instant  of  feeding. 


C  * 


Y(l)  =  ((C0Di-Y(l))*Q/2 
Y(2)  =  ((Vai-Y(2))*Q/2  + 
Y(3)  =  ((Xai-Y(3))*Q/2  + 
Y(4)  ■  ((Xmi-Y(4))*Q/2  + 
CH4T  =  0. 
C02T  =  0. 
GAST  =  0. 
4  05  CALL  RUNGE 


+  Y(l)*(V-Q/2))/V 
Y(2)*(V-Q/2))/V 
Y(3)*(V-Q/2))/V 
Y(4)*(V-Q/2))/V 


C  * 

C  ***  Update  time  and  CH4  production  in  this  printing  period 

C  * 

NT    =  NT  +  1 

TIME  =  NT  *  DELT  +  TINI 

TPD   =  TIME  +  DELT/2. 

CH4   =  V  *  Yd  *  MU2  *  Y(4)  /  Yc  *  DELT 

C02   =  (Ye*MUl*Y(3)  +  Yf*MU2*Y(4))  *  V  *  DELT 

GAS   =  CH4  +  CO  2 

CH4T  =  CH4T  +  CH4 

C02T  =  C02T  +  C02 

GAST  =  GAST  +  GAS 


C  * 
C  *** 
c  * 


404 


Time  to  output?  update  total  CH4  production 

IF  (AMOD(TPD,0.5)  .GE.  DELT)  GOTO  405 

IF  (AMOD(TPD,PRDEL)  .LE.  DELT)  GOTO  404 

Y(l)  =  ((CODi-Y(l))*Q/2  +  Y(l)*(V-Q/2))/V 

Y(2)  =  ((Vai-Y(2))*Q/2  +  Y(2) * (V-Q/2) )/V 

Y(3)  =  ((Xai-Y(3))*Q/2  +  Y(3) * (V-Q/2) )/V 

Y(4)  =  ((Xmi-Y(4))*Q/2  +  Y (4) * (V-Q/2) )/V 

GOTO  405 

SGPR  =  GAST/  ((RDAY-PDAY)*V) 

SMPR  =  CH4T  /  ((RDAY-PDAY)*V) 

PDAY  =  RDAY 

SPCH4  =  (CH4T  /  GAST)  *  100 

CALL  OUTPUT 


C 
C 
C 


***  Time  to  stop? 
* 

IF  (TPD  .LT.  TFIN)  GOTO  407 
CLOSE (46) 
CLOSE (47) 

CALL  SUMSQR (IBEGIN, IEND, SGPR,GPR,SSGPR) 
CALL  SUMSQR ( I BEGIN , IEND , SMPR , MPR , SSMPR) 
CALL  SUMSQR ( IBEGIN , IEND , SCOD , COD , SSCOD) 
CALL  SUMSQR ( IBEGIN , IEND , SPCH4 , PCH4 , SSPCH4 ) 
RSS  =  SSGPR  +  SSMPR  +  SSCOD  +  SSPCH4 


RETURN 
END 
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C  ********************************************************** 

C  *  * 

C  *                     #5  Init  * 

C  *  * 

C  *  Subroutine  to  initialize  values  of  state  variables  and  * 
C  *  assign  the  array  position  to  each  parameter. 

C  *  * 

C  ********************************************************** 

SUBROUTINE  INIT 

REAL*8  MU1 , MU2 , MUHAT1 , MUHAT2 , Ksl , Ks2 , Kdl , Kd2 , A ( 14 ) 
REAL  MP,MPR 

COMMON  /INPUT/NEQN,OUTFILE,NPARM 
COMMON  /Y0/COD0,VA0,XA0,XM0 

COMMON  /TIMER/ISAMP , TFIN , TINI , DELT , PRDEL , TIME , NT 
COMMON  /YVALUE/Ya,Yb,Yc,Yd,Ye,Yf 
COMMON  /KVALUE/Ksl,Ks2,Kdl,Kd2 
COMMON  / GROWTH/MUHAT 1 , MUHAT 2  ,  MU  1 ,  MU  2 
COMMON  /GAS/CH4T , GAST , MP , GP , GPR , MPR , PCH4 , 
&  SGPR,SMPR,SPCH4 

COMMON  /PARM/A,IC 


c  * 

c  *** 

Reset 

for  new  run 

c  * 

CH4T 

=  0. 

GAST 

=  0. 

GP 

=  0. 

MP 

=  0. 

NT 

=  0 

TIME 

=  TINI 

c  * 

c  *** 

Assign  array  position  to  each  parameter 

c  * 

Ya 

=  A(l) 

Yb 

=  A(2) 

Yc 

=  A(3) 

Yd 

=  A(4) 

Ye 

=  A(5) 

Yf 

=  A(6) 

Ksl 

=  A(7) 

Ks2 

=  A(8) 

Kdl 

=  A(9) 

Kd2 

=  A(10) 

MUHAT1  =  A (11) 
MUHAT2  =  A (12) 
XAO  =  A  (13) 
XMO     =  A  (14) 

RETURN 
END 
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C  ********************************************************** 

C  *  * 

C  *                     #6  Runge  * 

C  *  * 

C  *  Subroutine  to  use  4th  order  Runge-Kutta  method  for  * 

C  *  integration.  * 

C  *  * 

C  *    Y(t+dt)  -  Y(t)  +  dt  *  [Ydot(t)/6  +  PYdotl (t+dt/2)/3 ]  * 

C  *              +  PYdot2(t+dt/2)/3  +  PYdot3 (t+dt)/6]  * 

C  *  * 

C  *  where  * 

C  *    PYdotl (t+dt/2 ) ,  PYdot2(t+dt/2)  and  PYdot3 (t+dt)  are  * 

C  *   the  rate  found  by  using  PYl(t+dt/2),  PY2 (t+dt/2)  and  * 

C  *   PY3(t+dt)  respectively.  * 

C  *   and  * 

C  *    PY1 (t+dt/2)  =  Y(t)  +  dt/2  *  Ydot(t)  * 

C  *    PY2 (t+dt/2)  =  Y(t)  +  dt/w  *  PYdotl (t+dt/2)  * 

C  *    PY3(t+dt)    =  Y(t)  +  dt  *  PYdot 2 (t+dt/2)  * 

C  *  * 

C  *  Need  call  subroutine:  RATE  * 

C  *  * 
C  ********************************************************** 

SUBROUTINE  RUNGE 

DOUBLE  PRECISION  MU1 ,MU2 ,MUHAT1,MUHAT2 , Ksl , Ks2 ,Kdl, Kd2 
REAL  MP,MPR 

COMMON  /INPUT/NEQN , OUTFILE , NPARM 
COMMON  /Y0/COD0,VA0,XA0,XM0 
COMMON  /SVAR/Y(10) ,Ydot(10) ,COD,Va 
COMMON  /RUNKU1/PY1(10) ,PY2(10) ,PY3(10) 
COMMON  /RUNKU2/PYdotl(10) ,PYdot2(10) ,PYdot3(10) 
COMMON  /TIMER/ ISAMP , TFIN , TINI , DELT , PRDEL , TIME , NT 
COMMON  /YVALUE/Ya,Yb,Yc,Yd,Ye,Yf 
COMMON  /KVALUE/Ksl,Ks2,Kdl,Kd2 
COMMON  /GROWTH/MUHAT1 , MUHAT2 , MU1 , MU2 
COMMON  /TANK/Q,V,HRT,VIN 
COMMON  /INFLOW/ CODi,Vai,Xai,Xmi 
COMMON  /GAS/CH4T,GAST,MP,GP,GPR,MPR,PCH4, 
&  SGPR,SMPR,SPCH4 

C  * 

c  ***  Predict  PY1 

C  * 

CALL  RATE(Y,Ydot) 

DO  601  I=1,NEQN 

PY1(I)  =  Y(I)  +  DELT/2  *  Ydot(I) 
601  CONTINUE 

C  * 

C  ***  Predict  PYdotl  from  RATE 

C  * 

CALL  RATE (PY1, PYdotl) 
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C  * 

C  ***  Predict  PY2 

C  * 

DO  602  I=1,NEQN 

PY2(I)  =  Y(I)  +  DELT/2  *  PYdotl(I) 

602  CONTINUE 

C  * 

C  ***  Predict  PYdot2  from  RATE 

C  * 

CALL  RATE(PY2,PYdot2) 

C  * 

C  ***  Predict  PY3 

C  * 

DO  603  I=1,NEQN 

PY3(I)  =  Y(I)  +  DELT  *  PYdot2(I) 

603  CONTINUE 

C  * 

C  ***  Predict  PYdot3  from  RATE 

C  * 

CALL  RATE(PY3,PYdot3) 

C  * 

C  ***  Update  state  variables 

C  * 

DO  604  I=1,NEQN 

Y(I)  =  Y(I)  +  DELT  *  (Ydot(I)/6  +  PYdotl(I)/3 

&  +  PYdot2(I)/3  +  PYdot3(I)/6) 

604  CONTINUE 

RETURN 
END 


C  ********************************************************** 
c  *  * 

C  *  #7  Output  * 

C  *  * 

C  *   Subroutine  to  print  out  values  of  state  variables.    * 
C  *  * 

C  ********************************************************** 

SUBROUTINE  OUTPUT 

DOUBLE  PRECISION  MU1 ,MU2 ,MUHAT1 ,MUHAT2 , Ksl , Ks2 , Kdl, Kd2 

REAL  MP,MPR 

CHARACTER  OUTFILE*12 

COMMON  /INPUT/NEQN,OUTFILE,NPARM 

COMMON  /Y0/COD0,VA0,XA0,XM0 

COMMON  /SVAR/Y(10) ,Ydot(10) ,COD,Va 

COMMON  /TIMER/ ISAMP , TFIN , TINI , DELT , PRDEL , TIME , NT 

COMMON  /GAS/CH4T , GAST , MP , GP , GPR , MPR , PCH4 , 
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&  SGPR,SMPR,SPCH4 

WRITE (46, 2701) ISAMP, TIME, (Y(ISTAT) , ISTAT=1,NEQN) ,COD, 
&   VA,GP,MP,GPR,MPR,PCH4,GAST,CH4T,SGPR,SMPR,SPCH4 
2701  F0RMAT(1X,I3,1X,F6.2,6(1X,F6.2) , 2 (2 (IX, F7 . 2) , 
&  2(1X,F7.3) ,F7.2)) 

RETURN 
END 

C  ********************************************************** 
c  *  * 

C  *  #8  Sumsqr  * 

C  *  * 

C  *  Subroutine  to  find  the  residual  sum  square  between  the  * 
C  *  simulated  and  observed  data.  The  result  will  be  used  * 
C  *  to  evaluate  the  fitness  of  parameters.  * 

C  *  * 

C  ********************************************************** 

SUBROUTINE  SUMSQR (M, N, X, XX, S) 

REAL*8  R(16) ,G(16,14) ,A(14) ,MAXIM(14) ,MINIM(14) ,MX1, 
&   MX2 , MI , SUMNEW ( 14 ) , RSS , 
&   MUl,MU2,MUHATl,MUHAT2,Ksl,Ks2, 
&   Kdl,Kd2 

CHARACTER  OUTFILE*12 

REAL  MP,MPR 

COMMON  /SVAR/Y(10) ,YDOT(10) ,COD,VA 

COMMON  /TIMER/ISAMP , TFIN , TINI , DELT , PRDEL, TIME , NT , 

COMMON  /GAS/CH4T , GAST , MP , GP , GPR , MPR , PCH4 , 
&  SGPR,SMPR,SPCH4 

COMMON  /SS/SSGPR, SSMPR, SSCOD, SSVA, SSPCH4 ,RSS , SCOD, 
&  SVA,SXA,SXM 

OPEN (UNIT=4  6 , FILE= ■ CSIMUPE . OUT ' ) 

AX  =  0. 

DO  801  I=M,N 

READ(46,1801)  ISAMP, TIME, SCOD, SVA, SXA, SXM, 
&     COD,VA,GP,MP,GPR,MPR,PCH4,SGP,SMP, 
&     SGPR,SMPR,SPCH4 
1801   F0RMAT(1X,I3,1X,F6.2,6(1X,F6.2) , 2 (2 (IX, F7 . 2) , 
&         2(1X,F7.3) ,F7.2) ) 

AX  =  AX  +  ((X-XX)/XX)  **  2 
801  CONTINUE 
S  =  AX 
CLOSE (46) 

RETURN 
END 

C  ********************************************************** 
c  *  * 

C  *  #9  Rate  * 

c  *  * 

C  *       Subroutine  to  describe  the  model.  * 
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C  * 

C  *  User  is  allowed  to  have  up  to  10  rate  equations  for 

C  *  his  model. 

C 

c 
c 
c 
c 
c 
c 
c 


Notations: 

X(l)  =  Chemical  oxygen  demand,  g  COD/L  (COD) 

X(2)  =  Volatile  fatty  acid,  g  VA/L  (Va) 

X(3)  =  Acetogen,  g/L  (Xa) 

X(4)  =  Methanogen,  g/L  (Xm) 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 


********************************************************** 

SUBROUTINE  RATE(X,Xdot) 

DOUBLE  PRECISION  MU1 ,MU2 ,MUHAT1 ,MUHAT2 , Ksl,Ks2 ,Kdl ,Kd2 

DIMENSION  X(10) ,Xdot(10) 

COMMON  /YVALUE/Ya,Yb,Yc,Yd,Ye,Yf 

COMMON  / RVALUE/ Ksl,Ks2,Kdl,Kd2 

COMMON  /GR0WTH/MUHAT1,MUHAT2,MU1,MU2 

COMMON  /TANK/Q,V,HRT,VIN 

COMMON  / INFLOW/ CODi,Vai,Xai,Xmi 


MU1 
MU2 

Xdot(l) 
Xdot ( 2 ) 

Xdot ( 3 ) 
Xdot (4) 

RETURN 
END 


MUHAT1  *  X(l)  /  (Ksl  +  X(l)) 

MUHAT2  *  X(2)  /  (Ks2  +  X(2)) 

Q/V  *  (CODi  -  X(l))  -  MUl*X(3)/Ya 

Q/V  *  (Vai  -  X(2))  +  Yb*MUl*X(3)/Ya  - 

MU2*X(4)/Yc 

Q/V  *  (Xai  -  X(3))  +  MU1*X(3)  -  Kdl*X(3) 

Q/V  *  (Xmi  -  X(4))  +  MU2*X(4)  -  Kd2*X(4) 
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Optimizer 


********************************************************** 


* 
* 
* 
* 
* 
* 
* 


Filename:  OP. FOR 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


* 
* 
* 


An  interval  elimination  search  method  was  used  in  this  * 
program  to  find  the  value  of  the  optimal  control  * 
variable  for  a  nonlinear  self-tuning  regulator  * 
operated  at  a  set-point  output.   For  maximum  MPR       * 

C  *  operation,  objective  function  need  be  modified. 

C  *  Beginning  and  ending  time  need  be  defined  in  the  DATA 

C  *  statement:  TINI,  TFIN. 


Subroutine  (or  Function)  No, 

1.  INIT 

2.  FIBON 

3.  CALY 

4.  RUNGE 

5.  RATE 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 

Subroutines  RUNGE  and  RATE,  which  are  the  same  as  in   * 
program  PE.FOR,  will  not  be  listed  here.  * 

* 

Subroutines  called  in  main  program:  INIT,  FIBON        * 

* 

Substrate  used:  COD  * 

* 

Variables  (not  including  those  defined  in  Nomenclature  * 
and  PE.FOR) : 
ALPHA  =  Desired  accuracy 
FIB   =  Fibonacci  numbers 

QHIGH  =  Upper  bound  of  the  range  of  flow  rate  to 
search,  L/day 
=  Lower  bound  of  the  range  of  flow  rate  to 

search,  L/day 
=  Predicted  manipulated  flow  rate,  L/day 
=  Predicted  MPR  for  operating  at  QM,  L  CH4/L- 

day 
=  Set  point  of  MPR,  L  CH4/L-day 


QLOW 

QM 
YM 

Yst 


* 

* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 


********************************************************** 


c 
c 
c 


***  Main  Program 
* 

REAL*8  MU1 , MU2 , MUHAT1 , MUHAT2 , Ksl , Ks2 , Kdl , Kd2 , 
&        YM,Yst,SMPR 
REAL  MP,MPR,FIB(50) 
CHARACTER  0UTFILE*12 
COMMON  /INPUT/Y0(10) , NEQN, OUTFILE 
COMMON  /SVAR/Y(10) ,Ydot(10) ,COD,Va 
COMMON  /RUNKU1/PY1(10) ,PY2(10) ,PY3(10) 
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COMMON  /RUNKU2/PYdotl(10) ,PYdot2(10) ,PYdot3(10) 
COMMON  /TIMER/ISAMP , TFIN , TINI , DELT , PRDEL , TIME , NT 
COMMON  /YVALUE/Ya,Yb,Yc,Yd,Ye,Yf 
COMMON  /KVALUE/Ksl,Ks2,Kdl,Kd2 
COMMON  /GR0WTH/MUHAT1,MUHAT2,MU1,MU2 
COMMON  /TANK/Q,V,HRT,VIN 
COMMON  /INFLOW/CODi,Vai,Xai,Xmi 
COMMON  /GAS/CH4T,GAST,MP,GP,GPR,MPR,PCH4, 
&  SGPR,SMPR,SPCH4 

COMMON  /SET/Yst 
COMMON  /FIB/ALPHA,QLOW,QHIGH 
COMMON  /MANI/QM,YM 
DATA  V/10.0/ 

DATA  CODi/30 . 2595/ , Vai/0 . 37875/ , Xai/0 . 0/ , Xmi/0 . 0/ 
DATA  Y0(1)/2.751/,Y0(2)/0.120/ 
DATA  TINI/ 3 4 7. /, DELT/ 0.01/ 
DATA  ALPHA/. 001/ 

NEQN    =  4 
I COUNT  =  4 

C  * 

C  ***  Repeat  finding  local  optimum  at  different  starting 
C  ***  points  for  icount  times,  icount  can  be  changed  as 
C  ***  user's  wish.   The  range  of  flow  rate  to  search  is 
C  ***  different  for  each  iteration,  the  lower  bound  is  set 
C  ***  at  0  L/day  and  the  upper  bound  will  be  varied  from 
C  ***  0  to  1  L/day. 
C  * 

DO  1  11=1, ICOUNT 

PRINT  *,  'II  =  • ,11 

PRINT  *,  '    ' 

QLOW   =0.0 

QHIGH  =  1.0  -  (II-1)*0.25 

CALL  INIT 

CALL  FIBON 
1  CONTINUE 

STOP 
END 

C  ********************************************************** 
c  *  * 

C  *  #1  Init  * 

C  *  * 

C  *  Subroutine  to  input  the  parameters  and  initial         * 
C  *  conditions  of  state  variables.  * 

C  *  * 

C  ********************************************************** 


SUBROUTINE  INIT 

REAL*8  MU1 , MU2 , MUHAT1 , MUHAT2 , Ksl , Ks2 , Kdl , Kd2 , 
&        YM,Yst,SMPR 
REAL  MP,MPR,FIB(50) 
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CHARACTER  0UTFILE*12 

COMMON  /INPUT/Y0(10) ,NEQN, OUTFILE 

COMMON  /SVAR/Y(10) ,Ydot(10) ,COD,Va 

COMMON  /TIMER/ISAMP , TFIN , TINI , DELT , PRDEL , TIME , NT 

COMMON  /YVALUE/Ya,Yb,Yc,Yd,Ye,Yf 

COMMON  /KVALUE/Ksl,Ks2,Kdl,Kd2 

COMMON  /GROWTH/MUHAT1 , MUHAT2 , MU1 , MU2 

COMMON  /GAS/CH4T , GAST , MP , GP , GPR , MPR , PCH4 , 
&  SGPR,SMPR,SPCH4 

COMMON  /SET/Yst 

DATA  TFIN/419.5/,Yst/.40/ 
C  * 

C  ***  Read  parameters  found  in  PE.FOR 
C  * 

OPEN  (UNIT=31,FILE='PARM.OPT') 

READ  (31,*)  Ya,Yb,Yc,Yd,Ye,Yf ,  Ksl,Ks2  ,  Kdl,Kd2  , 
&  MUHAT1,MUHAT2,Y0(3) ,Y0(4) 

CLOSE (31) 

C  * 

C  ***  update  state  variables 

C  * 

DO  101  I=1,NEQN 

Y(I)  =  Y0(I) 
101  CONTINUE 


C 

* 

c 

*** 

Reset 

for  new  run 

C 

* 

CH4T 

= 

0. 

GAST 

= 

0. 

GP 

= 

0. 

MP 

= 

0. 

NT 

= 

0 

TIME 

= 

TINI 

RETURN 

END 

c  ********************************************************** 

c  *  * 

C  *                       #2  Fibon  * 

C  *  * 

C  *  Subroutine  to  optimize  the  control  variable  using  the  * 

C  *   Fibonacci  method.  * 

C  *  * 

C  *  Need  call  subroutine:  CALY  * 

C  *  * 
C  ********************************************************** 

SUBROUTINE  FIBON 

REAL*8  MU1 , MU2 , MUHAT1 , MUHAT2 , Ksl , Ks2 , Kdl , Kd2 , 
&        YM,Yst,SMPR 
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REAL  MP,MPR,FIB(50) 
CHARACTER  0UTFILE*12 
COMMON  /INPUT/Y0(10) ,NEQN,OUTFILE 
COMMON  /SVAR/Y(10) ,Ydot(10) ,COD,Va 
COMMON  /RUNKU1/PY1(10) ,PY2(10) ,PY3(10) 
COMMON  /RUNKU2/PYdot 1 ( 10 ) , PYdot2 ( 10 ) , PYdot 3 ( 10 ) 
COMMON  /TIMER/ISAMP, TFIN , TINI , DELT , PRDEL, TIME , NT , 
COMMON  /YVALUE/Ya,Yb,Yc,Yd,Ye,Yf 
COMMON  /KVALUE/Ksl,Ks2,Kdl,Kd2 
COMMON  /GROWTH/MUHAT1 , MUHAT2 , MU1 , MU2 
COMMON  /TANK/Q,V,HRT,VIN 
COMMON  /INFLOW/CODi,Vai,Xai,Xmi 
COMMON  /GAS/CH4T,GAST,MP,GP,GPR,MPR,PCH4, 
&  SGPR,SMPR,SPCH4 

COMMON  /SET/Yst 
COMMON  /FIB/ALPHA, QLOW, QHIGH 
COMMON  /MANI/QM,YM 

C  * 

C  ***  ALPHA  is  the  desired  accuracy,  QLOW  and  QHIGH  are  the 

C  ***  boundary  values  of  the  flow  rate  Q  and  DEL  is  the 

C  ***  range 

C  * 

DEL  =  QHIGH  -  QLOW 

C  * 

C  ***  Define  the  first  three  Fibonacci  numbers 

C  * 

FIBO    =1.0 

FIB(l)  =  1.0 

FIB(2)  =  2.0 

C  * 

C  ***  Calculate  the  remaining  Fibonacci  numbers 
C  * 

BB  =  1.0  /  ALPHA 

JJ  =  2 
3  01  JJ  -  JJ  +  1 

FIB(JJ)  =  INT(FIB(JJ-1) )  +  INT(FIB(JJ-2) ) 

CC  =  FIB(JJ) 

IF  (CC  .LT.  BB)  GOTO  3  01 


c 

* 

C 

*** 

First  step 

C 

* 

1  =  0 

KK  =  JJ  -  2 

IK  =  JJ  -  2 

BL  =  QHIGH  -  QLOW 

ALL  =  FIB(IK)*BL  /  FIB(JJ) 

Ql  =  QLOW  +  ALL 

Q2  =  QHIGH  -  ALL 

CALL  CALY(Q1,SMPR1) 

CALL  CALY(Q2,SMPR2) 
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JK  =  1 

C  * 

C  ***  Succeeding  steps 

C  * 

IK  =  IK  -  1 
JJ  =  JJ  -  1 
DO  302  1=1, KK 

IF  (SMPR2  .LE.  SMPR1)  THEN 
QLOW  =  QLOW  +  ALL 
BL  =  QHIGH  -  QLOW 
Ql  =  Q2 

CALL  CALY(Q1,SMPR1) 
ALL  =  FIB(IK)*BL  /  FIB(JJ) 
Q2  =  QHIGH  -  ALL 
CALL  CALY(Q2,SMPR2) 
II  =  I  +  1 
IK  =  IK  -  1 
JJ  =  JJ  -  1 
IF  (IK  .LT.  1)  IK=1 
ELSE 

QHIGH  =  QHIGH  -  ALL 
BL  -  QHIGH  -  QLOW 
Q2  =  Ql 

CALL  CALY(Q2,SMPR2) 
ALL  =  FIB(IK)*BL  /FIB(JJ) 
Ql  =  QLOW  +  ALL 
CALL  CALY(Q1,SMPR1) 
II  =  I  +  1 
IK  =  IK  -  1 
JJ  =  JJ  -  1 
IF  (IK  .LT.  1)  IK=1 
ENDIF 
302  CONTINUE 

C  * 

C  ***  Calculation  of  the  final  range  of  the  flow  rate 

C  * 

EPS  =  0.001  *  Ql 
DL   =  Ql  +  EPS 
CALL  CALY(DL,SMPR3) 
IF  (SMPR3  .LE.  SMPR2)  THEN 
CALL  CALY (QHIGH, SMPRH) 
QM  =  (Ql+QHIGH)/2 
ELSE 

CALL  CALY (QLOW, SMPRL) 
QM  =  (Ql+QLOW)/2 
ENDIF 

CALL  CALY(QM,YM) 
C  * 

C  ***  Check  COD  constraint,  Y(l)  <=  4.5  (85%  removal) 
C  * 

IREDUCE  =  0 
304  IF  (Y(l)  .LE.  4.5)  THEN 
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GOTO  303 
ELSE 

PRINT  *,  'COD  =  '  ,Y(1),'  g/L1 

QM  =  QM  -  0.05 

CALL  CALY(QM,YM) 

GOTO  304 
ENDIF 

C  * 

C  ***  Print  the  results 

C  * 

303  WRITE(6,2307)  QM,SMPR,Y(1) 
2307  FORMAT (/ IX, 'Manipulated  flow  rate  QM  =  ',F10.4,' 
&  L/day'/ IX, 'Predicted  output  YM      =  ',F10.4,' 
&  L/L/day'/lX, 'Predicted  COD  =  ',F10.4) 

RETURN 
END 

C  ********************************************************** 

c  *  * 

C  *  #3  Caly  * 

C  *  * 

C  *  Subroutine  to  conduct  the  dynamic  simulation  of  * 
C  *  methane  fermentation  and  find  the  residual  square  of  * 
C  *  the  desired  output  ([Yest  -  Yst]  *  2)  as  the  objective  * 
C  *  function  used  in  subroutine  FIBON.  The  fourth  order  * 
C  *  Runge-Kutta  method  is  applied  here  for  integration.  * 
C  *  * 

C  *  Need  call  subroutines:  RUNGE  * 

C  *  * 

C  ********************************************************** 

SUBROUTINE  CALY(QF,SY) 

REAL*8  MU1 , MU2 , MUHAT1 , MUHAT2 , Ksl , Ks2 , Kdl , Kd2 , 
&        YM,Yst,SMPR 

REAL  MP,MPR 

COMMON  /INPUT/Y0(10) ,NEQN,OUTFILE 

COMMON  /SVAR/Y(10) ,Ydot(10) ,COD,Va 

COMMON  /RUNKU1/PY1(10) ,PY2(10) ,PY3(10) 

COMMON  /RUNKU2/PYdotl(10) ,PYdot2(10) ,PYdot3(10) 

COMMON  /TIMER/ISAMP,TFIN,TINI,DELT,PRDEL,TIME,NT 

COMMON  /YVALUE/Ya,Yb,Yc,Yd,Ye,Yf 

COMMON  /KVALUE/Ksl,Ks2,Kdl,Kd2 

COMMON  /GROWTH/MUHAT1 , MUHAT2 , MU1 , MU2 

COMMON  /TANK/Q,V,HRT,VIN 

COMMON  / INFLOW/ CODi, Vai , Xai, Xmi 

COMMON  /GAS/CH4T,GAST,MP,GP,GPR,MPR,PCH4, 
&  SGPR,SMPR,SPCH4 

COMMON  /SET/Yst 

PQ     =  0.452 

NT     =0 

TPD    =  TINI  +  DELT/2 


OPEN  (UNIT=37,FILE='COPT.DATl) 
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C 
C 
C 
C 
C 
C 
C 
C 
C 


***  Given  N  day's  data,  CALY  is  able  to  predict  N+0.5 

***  day's  data  except  flow  rate  Q.   By  using  the  QF  from 

***  FIBON,  we  may  predict  the  N+l  day's  output. 

***  Therefore,  a  0.5  day  mark  is  set  in  next  line. 

***  This  program  is  designed  to  predict  the  next 

***  24  hours  performances  using  the  predicted  QM  for  next 

***  12  hours. 

407  IF  (TPD  .GT.  TFIN)  THEN 

Q      =  PQ 
ISAMP  =  ISAMP  +  1 
RDAY   =  RDAY  +0.5 
IF  (TPD  .GT.  (TFIN+0.5))  Q  =  QF 
GOTO  4  02 
ENDIF 

Q  -  PQ 

READ(37,1401)  ISAMP, RDAY, PQ,VIN,GP, MP, GPR,MPR, PCH4 , 

&  COD,VA 

1401  FORMAT(1X,I3,F6.1,2F6.3,2F6.2,2F6.3,F6.2,2F7.0) 

COD  =  COD  /  1000. 
VA    =  VA  /  1000. 


C  * 

C  ***  Finding  the  concentration  of  state  variables  at  the 

C  ***  instant  of  feeding. 

C  * 


402  Y(l)  =  ((CODi-Y(l))*Q/2 
Y(2)  =  ((Vai-Y(2))*Q/2  + 
Y(3)  =  ( (Xai-Y(3) ) *Q/2  + 
Y(4)  =  ((Xmi-Y(4))*Q/2  + 

403  CH4T  =  0. 
C02T  =  0. 
GAST   =  0. 

4  05  CALL  RUNGE 


+  Y(l)*(V-Q/2))/V 
Y(2)*(V-Q/2))/V 
Y(3)*(V-Q/2))/V 
Y(4)*(V-Q/2))/V 


C  * 

C  ***  Update  time  and  CH4  production  in  this  printing  period 

C  * 

NT    =  NT  +  1 

TIME  =  NT  *  DELT  +  TINI 

TPD   =  TIME  +  DELT/2. 

CH4   =  V  *  Yd  *  MU2  *  Y(4)  /  Yc  *  DELT 

C02   =  (Ye*MUl*Y(3)  +  Yf*MU2*Y(4))  *  V  *  DELT 

GAS   =  CH4  +  C02 

CH4T  =  CH4T  +  CH4 

C02T  =  C02T  +  C02 

GAST  =  GAST  +  GAS 
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C 
C 
C 


* 

*** 

* 


Time  to  print?  update  total  CH4  production 


404 


IF 

IF 

Y(l) 

Y(2) 

Y(3) 

Y(4) 

GOTO 

SGPR 

SMPR 

SPCH4 


(AMOD(TPD,0.5) 
(AMOD(TPD,RDAY) 


GE.  DELT)  GOTO  405 
.LE.  DELT)  GOTO  4  04 


=  ((CODi-Y(l))*Q/2  +  Y(l)*(V-Q/2) )/V 
=  ((Vai-Y(2))*Q/2  +  Y (2) * (V-Q/2) )/V 


=  ((Xai-Y(3))*Q/2  + 
=  ((Xmi-Y(4))*Q/2  + 
405 
=  GAST  /  (0.5*V) 
=  CH4T  /  (0.5*V) 
=  (CH4T  /  GAST)  *  100 


Y(3)*(V-Q/2))/V 
Y(4)*(V-Q/2))/V 


C  * 

C  ***  Time  to  stop? 

C  * 

IF  (TPD  .GT.  (TFIN+1.0))  GOTO  406 
IF  (TPD  .LT.  (TFIN+1.0))  GOTO  407 
406  CLOSE(37) 

C  * 

C  ***  Find  value  of  the  objective  function 

C  * 

SY  -  (SMPR  -  Yst)  **  2 

RETURN 
END 
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Control  Program 


10 

20 

30 

40 

50 

60 

70 

72 

74 

76 

78 

80 

82 

84 

86 

400  COLOR  7,1:CLS 

410  ' 

500  'input  the  flow  rate  Q  and  calculate  the  pumping  time 

510  'pumping  time  for  reactor  C  =  cpmin  minutes  and  cpsec 

seconds. 
520  • 
530  INPUT  "the  optimal  flow  rate  for  reactor  C  (in 

milliliter)  =  ",CQ 
540  • 
550  'volume  to  feed  for  next  12  hours  will  be  egual  to  Q/2 

Liter 
560  ' 

570  'pumps  were  operated  at  a  constant  speed  of  15  ml/min 
580  ' 

590  CPMIN=INT(CQ/(2*15) ) 

600  CPSEC=INT((CQ/ (2*15) -CPMIN) *60+0. 5) 
700  ' 
1000  CLS 
1002  LOCATE  1,15:PRINT  "pumping  time  for  reactor  C  is 

"; CPMIN;"  minutes" ; CPSEC ; "  seconds" 
1006  LOCATE  16, 25: COLOR  14: PRINT  "Program  is  running 


and  binary  code  for  relay  #1 


Program  Name:  CONTROL. BAS 

This  program  was  designed  to  use  relay  to  control  the 
power  of  influent  and  effluent  pumps,  and  calculate  the 
pumping  time  from  the  input  flow  data  found  in  the 
optimizer. 

relay  #1:  binary  code  =  2  for  power  on. 

INIT,  IONAME  and  DIGWRITE  are  functions  defined  in 
SOFT500  which  is  a  software  specially  designed  for 
Keithley  System  500  series. 


1007 

COLOR  7,1,1 

1008 

■ 

1009 

1  reset  swit 

1010 

i 

1012 

CCODE=0 

1020 

SWITCH=0 

1021 

i 

1025 

' initialize 

1026 

i 

1030 

CALL  INIT 

1031 

i 

1032 

'Set  up  the 

"7"  is  the 

control . 

system 


IONAME  command.   "relay"  is  the  task  nai 
System  570  slot  designation  for  power 
ra"  is  the  port  location  for  channel  0-7 
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1033  ' 

1040  CALL  IONAME' ("relay" , 7 , "a") 

1050  KEY  OFF 

1055  ' 

1060  'TIME$  =  HH:MM:SS 

1070  'HH$=LEFT$(TIME$,2) :MM$=LEFT$ (MID$ (TIME$ , 4) ,2) : 

SS$=RIGHT$ (TIME$ , 2 ) 
1080  ' 

1100  'check  if  it  is  time  for  feeding  the  reactor 
1110  IF  (LEFT$(TIME$,2)="09"  OR  LEFT$ (TIME$ , 2 ) ="21" )  AND 

LEFT$(MID$(TIME$,4) ,2)="05"  AND  RIGHT$ (TIME$ , 2) ="00" 

THEN  GOSUB  2000 
1120  ' 

1200  'check  if  it  is  time  to  stop  feeding  reactor  C 
1210  ' 
1220  IF  (LEFT$(TIME$,2)="09"  OR  LEFT$ (TIME$ , 2) ="21")  AND 

VAL(LEFT$(MID$(TIME$,4) ,2))=5+CPMIN  AND 

VAL(RIGHT$(TIME$,2) )=CPSEC  THEN  GOSUB  3  000 
1240  GOTO  1100 
1250  END 
1260  ' 

2000  •****  feeding  pumps  are  On  **** 
2005  ' 
2010  COLOR  4, 7: LOCATE  6 , 5 : PRINT  TIME$: LOCATE  6, 15: PRINT 

"start  feeding  for  reactor  C  ..." 
2020  CCODE=2 
2030  SWITCH=CCODE 
2040  ' 

2  050  'Do  the  digital  output.   "relay"  is  the  I/O  task  name. 

"switch"  is  the  output  value  which  controls  the  ON/OFF 

status  of  the  channels  on  port  A. 
2060  • 

2080  CALL  DIGWRITE' ("relay", switch) 
2090  RETURN 
2100  ' 

3000  •****  stop  feeding  for  reactor  C  *** 
3010  ' 
3020  COLOR  7,8:LOCATE  6 , 5 : PRINT  TIME$: LOCATE  6, 15: PRINT 

"stop  feeding  for  reactor  C  ...." 
303  0  CCODE=0 
3040  SWITCH=CCODE 
3050  ' 

3  060  'Do  the  digital  output. 
3070  ' 

3080  CALL  DIGWRITE' ("relay", switch) 
3  090  RETURN 
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